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Synopsis 


SYNOPSIS 

In this thesis, the effect of free surfaces on fracture characteristics for elastostatic < 
elastodynamic problems is investigated using three-dimensional boundary element meth 
Free surface effect is studied for a through center cracked finite body for three situatic 
quasi-static conditions, stationary crack under impact loading and running crack un 
impact loading. 

A crack front intersecting the free surface in a 3-D body has been a subject of intei 
for last three decades. Through-thickness fatigue cracks in plate geometries almost alw 
show crack front tunnelling, i.e., a curved crack front with the deepest point of the cr 
at the center. Benthem [1] analyzed a problem of a quarter-plane crack in a semi-infii 
half-space subjected to symmetric quasi-static loading. For a crack front normal to the 
surface, he derived that at the corner point, stresses have r ~ A 52 singularity for Poisst 
ratio (u) of 0.3. By a local finite element analysis, Bazant and Estenssoro [2] estimated 
critical intersection angle to be 101° (for v = 0.3) for a stationary crack under elastost 
loading. By a finite element formulation similar to that used by Bazant and Estenssoro 
Gudmundson and Ostulund [3] found that a weaker singularity also exists at the coi 
point for a dynamically growing crack and estimated the critical angle to be 101°, aln 
independent of the crack-tip velocity. They found that for angles between 90° and afc 
101°, singularity weakened further with increasing velocity and for angles larger than at 
101°, singularity became stronger with increasing velocity. 

The success enjoyed by these local numerical analyses is not found to the same exter 
finite element treatments of more global or complete problems which contain crack/ sur: 
intersections. Different predictions were made concerning the variation in energy ret 
rate ( G ) with some results showing an increase in G in the free surface region, son 
decrease and still others no change for quasi-static analysis. Burton et al. [4] generj 
critical crack profiles using FEM with each profile having less than 2% numerical variai 
in G with depth. The crack profiles generated by them were significantly less curved t 
those observed in practice [5]. They surmised that underlying assumtions in linear els 
fracture mechanics (LEFM) were to a degree incapable of explaining all of the var; 
phenomena leading to free surface retardation. Using FEM for analysis, Lin and Smitl 
found that slightly non-orthogonal intersection caused a great loss of path-independi 
in the value of J integral up to 47% near the free surface. Thus in global situation, fi 
element analysis is not able to obtain an accurate stress state adjacent to the crack 
in the region close to the free surface, especially for cracks which are not normal tci 
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surface. 

In this thesis, BEM is used to study the effect, of free surface on the state of stress 
at the crack vertex for a center cracked finite specimen under both static and dynamic 
conditions. Straight as well as curved crack fronts are analysed. Influence of Poisson’s 
ratio and intersection angle of crack front with free surface on stress intensity factor along 
the crack front is studied. 

Quasi-static analysis 

Earlier studies have concentrated on calculation of stress singularity exponent on point of 
intersection of crack front with free surface and decay in value of energy release rate near 
free surface. This study capitalizes on the strength of BEM to calculate J\ l integral more 
accurately near free surface as stresses and their gradients can be calculated more accurately 
using global boundary integral equations. Stress intensity factors can be calculated in BEM 
using traction values at the nodes along the crack front. A dose correspondence between 
Ki computed from traction values and Jj n integral is seen in this study. Thus, it is possible 
to relate the evaluated state of stress at the crack vertex, where crack front, meets the free 
surface, with a value of SIF. At free surface, reduction in value of ./" up to a value of 50% 
of plane-strain value is estimated. Thickness of boundary layer where free surface effects 
are predominant is estimated. Critical intersection angle of crack front with free surface 
computed in this work for a finite center cracked specimen matches closely with the value 
predicted by the local FEM analysis of Bazant, and Estenssoro for a half-space problem [2]. 

Stationary crack under impact loading 

The three-dimensional time-domain displacement boundary integral formulation is imple- 
mented for analysis of a center cracked finite solid under symmetrically applied tensile step 
loading. A new partitioning scheme is proposed for spatial integration of elastodynamic 
kernels. Algorithm for calculating principal value integrals suggested by Guiggiani and Gi- 
gante [7] is extended to elastodynamic analysis. Also, combination of projection technique 
with semi-analytical integration, used for nearly singular integrals for static problems by 
Cruse and Aithal [8] and Mi and Aliabadi [9], is extended to elastodynamic problems. 

In this thesis, variation of dynamic stress intensity factor (DSIF) along the crack front 
is studied. The state of stress is evaluated at the crack vertex, where crack front, meets 
the free surface. The effect of free surface on DSIF is investigated. The effect of waves 
travelling in thickness direction is explained. It is possible to estimate accurately the 
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critical intersection angle of the crack front with the free surface at which square-roc 
singularity is restored at the crack vertex under step loading. 

Running crack under impact loading 

The three-dimensional time-domain displacement boundary integral formulation is use 
for analysis of a moving center-crack in a finite solid under suddenly applied unifon 
normal pressure to crack faces. The algorithm used by Gallego and Dominguez [10] for 
running crack in 2-D analysis is extended to 3-D analysis in this thesis. The crack growt 
is modelled by using moving singular elements and a remeshing technique. For movir 
elements, space- and time-dependent interpolation functions are used. 

Variation of dynamic stress intensity factor (DSIF) along the crack front is studied. T1 
state of stress is evaluated at the crack vertex, where crack front meets the free surfac 
The effect of free surface on DSIF is investigated. It is again possible to estimate accurate 
the critical intersection angle of the crack front with the free surface at which square-ro< 
singularity is restored at the crack vertex under step loading for various speeds of cra< 
propagation. 

Main findings 

• In earlier studies by Benthem [1] and Bazant and Estenssoro [2], local eigenvali 
analysis was done for a quarter-plane crack in a semi-infinite half-space subject! 
to symmetric loading. This is the first study for a complete 3-D finite geomet 
problem that obtained the same critical intersection angle as obtained by Baza 
and Estenssoro [2] for the half-space problem. This thesis establishes that LEF 
can explain 101° angle as critical angle of intersection for a Mode-I finite specime] 

• Main highlight of the thesis is bringing out the relevance of SIF directly calculat 
from crack front tractions and its correspondence with J[ L integral values, especia] 
near the free surface. This approach may facilitate study of free surface effect i 
more different problems. In FEM analysis, there is loss of path-independence 
evaluation of J integral near the free surface. In this study, it is found that desp: 
some mesh dependence, SIF value at the crack vertex, calculated from traction 
the vertex, reflects the level of singularity at the vertex. Its value is sensitive 
changes in vertex singularity brought about by changes in Poisson’s ratio or angle 
intersection of crack front with the free surface for the same level of mesh refineme: 
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Chapter 1 
Introduction 


In the theory of elasticity, two types of two-dimensional problems are recognised: one of 
plane stress and one of plane strain. A state of plane stress is approached in a very thin bar 
without discontinuities in tension or compression. A state of plane strain is approached in 
the central transverse section of an infinitely long cylinder subjected to uniform loading 
along a line parallel to its axis. In dynamic situations, a plane strain case is obtained if 
the wavelength is small compared with all dimensions of the body of material in which the 
waves travel. A plane stress case is obtained if the waves are propagated in a plate such 
that the wavelength is large compared with the thickness of the plate. 

When a thin bar in tension has a discontinuity, such as a hole, groove or crack, the 
state of stress may be affected by the discontinuity. It is both of theoretical and practical 
interest to examine the effect of the thickness on the state of stress resulting from such 
discontinuities. Approximate theories of thin plates are unreliable in the case of plates of 
considerable thickness, or when the plates are weakened by cavities whose dimensions are 
of the same order of magnitude as the plate thickness. In such cases, the stress variations in 
the thickness direction must be accounted for and the problem adopts a three-dimensional 
character. 

The stress distribution in a thick plate containing a smooth circular cavity has been 
discussed by Sternberg and Sadowsky [13]. Inner layers are seen to carry a somewhat 
greater stress than the outer layers of the bar. Their work showed that the thickness of 
the plate can exert appreciable influence on the stress concentrations of the circular hole. 
When the periphery of the cavity contains re-entrant corners or singular points, such as 
a sharp crack, the problem is considerably more difficult mainly because the conventional 
mathematical techniques are not suitable for handling 3-D problems with geometric sin- 
gularities. Also, the stress distribution in a cracked plate of finite thickness is neither in a 
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state of plane stress or plane strain nor generalized plane stress. 

Hartranft and Sih [14] considered a through-thc- thickness crack into sec ting normally 
to the free surfaces in a finite size plate. They defined a boundary thickness /I near the free 
surface by an empirical equation 0 = t /( 4 + 16 t/a) in which the square-root singularity is 
not valid. In this equation, t is the thickness of the plate and 2a is the crack length. But 
they did not analyse the state of stress in boundary layer region near the crack vertex. 
Crack vertex is the corner point where the crack front intersects the free surface for a 
through-thickness crack. 

Yang and Freund [15], motivated by concerns regarding experimental measurements 
made in the vicinity of the crack tip, explored the 3-D crack problem using a boundary 
layer approach. They studied the state of stress in a thin plate which contained through 
cracks with a view toward assessing the influence of transverse shear on the crack tip stress 
and deformation fields. Their boundary layer solution showed a smooth transition from 
the near-tip generalized plane strain solution to plane-stress solution at distances from the 
crack tip V’ of 1/2 to 3/4 of plate thickness ‘/i’ depending on the value of Poisson’s ratio. 
There is a zone around the crack edge where three-dimensional effects influence the state 
of stress. They stated that the plane stress assumption is typically valid only for r > h/2. 

Rosakis and Ravi-Chandar [16] studied the extent of the region of three-dimensionality 
of the crack-tip field using reflected and transmitted caustics. They found that plane stress 
conditions prevail at distances from the crack-tip greater than half the specimen thickness, 
while no significant plane strain region was detected by them. 

Yang and Freund [15] and Rosakis and Ravi-Chandar [16] were concerned about the 
extent of zone of influence of three-dimensional effects around the crack edge so as to 
identify the regions in which local experimental measurements based on two-dimensional 
theory can be performed with confidence. But they did not study the effect, of free sur- 
face on state of stress at the crack vertex. When a crack front intersects a free surface, 
elastic analysis and experiment suggest that the inverse-square root singularity of LEFM 
is altered. Benthem [1] analyzed a problem of a quarter-plane crack in a semi-infinite half- 
space subjected to symmetric loading. He derived that at the corner point, stresses have 
t singularity for Poisson’s ratio v equal to 0.3. Although it does not seem a dramatic 
change, it implies that the stress intensity factor at this point is equal to zero (or does 
not exist) since it is defined as the strength of the square-root singularity at the crack tip. 
It also causes the energy release rate to be equal to zero, since this can have a bounded, 
non-zero value only in the case of an exact r~ l singularity of the strain energy. Bazant 
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and Estenssoro [2] confirmed Benthem’s values for the singular eigenvalue at the crack 
vertex through a 3-D finite element analysis developed for conical regions. In addition, 
they considered intersections other than normal and found angles at which square-root 
singularity is restored. They used extrapolation procedures to estimate singularity indices 
corresponding to infinite degrees of freedom from results of eigenvalue analysis of a set of 
uniform grids. By a finite element formulation similar to that used by Bazant and Es- 
tenssoro [2], Gudmundson and Ostulund [3] found that a weaker singularity also exists at 
the corner point for a dynamically growing crack and estimated the critical angle to be 101° 
for v = 0.3, almost independent of the crack-tip velocity. The success enjoyed by these 
local numerical analyses for half-space problems is not found to the same extent in finite 
element treatments of more global or complete problems (i.e., finite bodies) which contain 
crack/surface intersections. Different predictions were made concerning the variation in 
energy release rate ( G ) with some results showing an increase in G in the free surface 
region, some a decrease and still others no change. Burton et al. [4] used a special set of 
integral equations for an elastic half-space weakened by a crack of finite width and infinite 
depth subjected to transverse excitation. They confirmed the existence of a weaker sin- 
gularity at global crack/surface intersections with both integral equation formulation and 
finite element analysis of a finite geometry with adjustable singularity index at the crack 
vertex. They found by FEM a decay of G as free surface is approached [zja — ¥ 0) with the 
drop between the maximum value and that of the closest Gauss point (at z/a = 0.015) to 
the free surface being 30% of plane-strain value G v when Poisson’s ratio is 0.4. Nakamura 
and Parks [17] estimated thickness of boundary layer region, where free surface effects are 
predominant, as 3% of plate thickness. They used domain integral method to calculate 
local energy release rate J and calculated a drop in value of J of the order of 30% of two- 
dimensional plane-stress value of J near the crack vertex. In FEM, calculation of J-integral 
at various positions along the crack front is constrained by the fact that accurate values of 
stresses are available only at Gauss points of the elements. Hence accurate determination 
of thickness of boundary layer at free surface becomes difficult as it requires having very 
thin elements near the free surface. 

It can be observed in experiments that an initially straight-fronted crack in specimens 
such as a compact-tension specimen usually grows to a slightly curved shape under a fatigue 
load [5]. Using local eigenvalue analysis with FEM, Bazant and Estenssoro [2] found the 
critical intersection angle to be 101°, a value also observed experimentally [5] for v — .3. 
Three-dimensional LEFM implies that for v ^ 0 the straight profiles associated with the 
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varying energy release rates would not propagate as such but rather take up curved profiles 
with constant energy release rates. Burton et al. [4] generated critical crack profiles using 
FEM with each profile having less than 2% numerical variation in G with depth. The crack 
profiles generated by them were significantly less curved than those observed in practice [5], 
Similarly, using FEM, Bakker [18] predicted a relative tunnelling depth of 2.5 - 3% based 
on constant stress intensity distributions for self-similar crack propagation. Experimentally 
observed tunnelling depths on fatigue cracks are generally larger, closer to 5%, which is 
close to circular front. Using FEM for analysis, Lin and Smith [6] found that slightly non- 
orthogonal intersection caused a great loss of path-independence in the value of J-integral 
up to 47% near the free surface. Thus, in finite geometry situation, finite clement analysis 
is not able to obtain an accurate stress state adjacent to the crack front in the region close 
to the free surface, especially for cracks which are not normal to the surface. 

Burton et al. [4] surmised that the underlying theoretical assumptions and developments 
in LEFM were to a degree incapable of explaining all of the various phenomena actually 
present in retardation of crack front near the free surface. Bakker [18] mentioned the need 
of invoking plasticity and crack closure effects to explain free surface retardation. On the 
other hand, for quasi-brittle behaviour with only small-scale yielding, Smith [1.9] found 
that cracks identical in shape to those in fatigued metal models could be grown under 
monotonic loads ( R = ratio of minimum Kj to maximum Ki — 1.0) in photoelastic bodies 
suggesting that closure effects might not strongly affect crack shapes. 

The calculation of SIF in two dimensions (for virtually arbitrary configurations) is com- 
monplace today. In three-dimensions, the case is slightly different. Existing methodology 
is sufficient for the solution of Mode-I problems provided the intersection of crack fronts 
with free surfaces is avoided and provided the crack front does not exhibit abrupt changes 
in curvature. The problems of brittle fracture involving intersection of the crack fronts 
with the free surfaces still require extremely careful examination. 

The boundary element method (BEM) is already established as an economical and 
accurate numerical technique for elastic stress analysis. It is particularly suitable for prob- 
lems that are 3-D and involve high local stress gradients, such as those encountered near 
cracks. It is possible to obtain higher resolution of interior stresses and their derivatives 
with BEM. The high accuracy of BEM makes it possible to analyse accurately the effect 

of free surface on the state of stress at the point of intersection of crack front with the free 
surface in this work. 
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1.1 Overview of Linear Elastic Fracture Mechanics 

A number of books and monographs are available to provide background on fracture me- 
chanics from various perspectives. For example, the books by Broek [20], Kanninen and 
Popelar [21], Cherepanov [22], Hellan [23] and Kumar [24] develop the topic from the me- 
chanics perspective. The books by Knott [25] and Lawn [26] introduce the subject from the 
mechanics perspective, whereas Barsom and Rolfe [27] discuss fracture within the context 
of materials selection and structural applications. 

Fracture mechanics is concerned with the quantitative description of the mechanical 
state of a deformable body containing a crack or cracks with a view toward characterizing 
and measuring the resistance of materials to crack growth. 

There being multiple flaws in real material, the growth of the most dangerous of which 
leads to fracture. Griffith [28] explained why the strength of real materials is from 2 to 3 
orders of magnitude lower than the theoretical strength implied by the perfect crystal struc- 
ture of bodies. Though Griffith’s analysis referred to sub-molecular size of micro-cracks, 
his technique could be simply extended to any macroscopic level by suitable redefinition 
of specific constants like surface energy of the material. Irwin [29] and Orowan [30] ex- 
tended Griffith’s theory to micro-cracks in steel and proposed the concept of quasi-brittle 
fracture. They included the region of irreversible strains near a crack edge into the surface 
of fracture and introduced the concept of effective surface energy which includes plastic 
work. This made it possible to consider fractures within the framework of the model of an 
ideally brittle body. 

Continuum field approach to fracture of solids was launched with the introduction of the 
elastic stress intensity factor, K , as a crack tip field characterizing parameter by Irwin [31]. 
This idea provided an alternate framework for discussing the strength of cracked solids of 
nominally elastic material. Irwin proposed that a crack will begin to grow in a cracked 
body with limited plastic deformation when K is increased to a value called the fracture 
toughness of the material. 

The crack tip region is said to be autonomous at fracture initiation or during crack 
growth if 

1. The extent of the region of nonlinearity from the crack edge is very small compared 
to all other length dimensions of the body and loading system. 

2. The mechanical state within this end region at incipient growth or during growth is 
independent of loading and geometrical configuration. 
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For the particular case of an elastic-plastic material, the property of autonomy implies 
that the crack tip plastic zone is completely surrounded by an elastic SIF held and that 
the state of stress within the plastic zone is determined by the level ot stress intensity ol 
the surrounding field. This situation, which was termed small-scale yielding by Rice [32], 
is the situation in which SIF is a useful fracture characterizing parameter. 

It may be noted that a stress distribution with an infinite singularity is clearly a math- 


ematical idealization, in that no real material can actually support such a stress. The 
rationalization for admitting the singular stress distribution, the strength of which is mea- 
sured by the SIF, is based on the concept of small-scale yielding. It is thus assumed that in 
the immediate vicinity of the crack tip, the potentially large stresses are relieved by some 
nonlinear process in a region whose dimensions are small compared to crack length and 
body dimensions. It is assumed further that the stress distribution in the elastic material 
adjacent to the small zone is adequately described by the dominant singular term in the 
elasticity solution. Under small-scale yielding conditions, the SIF may In' considered to 
be a one-parameter measure of the amplitude of the stress which is being applied to the 
material in the crack tip region. The SIF approach circumvents consideration of how the 
material in the crack tip region actually responds to the applied stress. 

Fracture mechanics is devoted to determining SIF or energy release rate for various 
conditions. For Mode-I loading, the components of stresses near the crack edge can be 
represented asymptotically as 




Kj 

y/2- 


nr 


■■faW 


( 1 . 1 ) 


Asymptotic solutions show 


1. universal spatial dependence 

2. All information about loading and configuration are embedded in the scalar multiplier 
called SIF. 


1.2 Overview of Dynamic Fracture Mechanics 

A few comprehensive books and monographs on dynamic fracture mechanics, notably, are 
Freund [12], Parton and Boriskovsky [33,34] and Broberg [35]. Dynamic fracture mechanics 
is the subfield of fracture mechanics concerned with fracture phenomena for which the role 
of material inertia becomes significant. Inertial effects can arise either from rapidly applied 
loading on a cracked solid or from rapid crack propagation. In the case of rapid loading, 
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the influence of the loads is transferred to the crack by means of stress waves through th< 
material. To determine whether or not a crack will advance due to the stress wave loading 
it is necessary to determine the transient driving force on the crack. 

Quasi-static, fracture mechanics gives the answer to the question of whether or not th< 
crack will grow catastrophically. The scope of dynamic fracture mechanics is much wider 
Quasi-static fracture mechanics formulates only the criterion of unstable crack propagation 
whereas dynamic fracture mechanics requires the formulation of a series of criteria — those 
for crack initiation, crack propagation, crack arrest, bending of crack trajectory, and cracl 
branching. Some basic problems of dynamic fracture mechanics can be listed as follows: 

1. Establishment of the frequency dependence of stress intensity factor for a crack sub- 
jected to harmonic loading. 

2. Establishment of temporal dependence of stress intensity factors for a crack undei 
impact loading. 

3. Establishment of the dependence of stress intensity factors on time and crack prop- 
agation velocity. 

4. Establishment of the law of crack propagation for a known dependence of the stress 
intensity factors on the crack propagation velocity. 

Theoretical treatments by Yoffe, Broberg and others show that the nature of the near 
field about the tip of a propagating crack changes with velocity (see Freund [12]). The 
stresses can be represented in the form analogous to formula (1.1) with the only difference- 
now the stress intensity factors depend on time whereas the angular stress distribution 
depends on the velocity, i.e., 

( 1 . 2 ) 

V 27rr 

The stress intensity factor in an impact loading substantially exceeds its static value 
in the transient period. Freund, in his series of papers [11,36], established the univer- 
sal analytic methods for evaluation of the transient stress intensity factor of a crack in a 
two-dimensional geometrical configuration under dynamic loading. One of his major con- 
clusions was that the stress intensity factor of a half-plane running crack is given by the 
product of the universal function of crack-tip velocity and the stress-intensity factor for an 
equivalent stationary crack [36]. 

The most popular model of dynamic fracture mechanics deals with growth of a rectilin- 
ear crack in an elastic plane. It is assumed that the energy consumed during the formation 
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of the unit area of a new surface is a material characteristic. The elastodynamic stress 
field at the crack tip is calculated to formulate the criterion of crack propagation as an 
equation of energy balance. 

Thus, the idealized model of dynamical fracture mechanics is based on two main prin- 
ciples: 

1. Stress fields at the crack tip are described in terms of stress intensity factors. 

2. The criteria of the crack start, arrest and propagation are derived from the energy 
balance equation. 

Thus, to verify the accuracy of the model to a real process, it is necessary to check the 
validity of these two statements. With this aim, one can analyse the experimental data on 
stress intensity factors and then compare the necessary experimental conditions for crack 
initiation, propagation and arrest with the corresponding theoretical predictions. 

In most experiments, correspondence between the SIFs and the crack propagat ion ve- 
locity turns out to be non-unique. At low loading rates and moderate loads, there is 
good agreement between the theoretical and experimental stress intensity factors, which is 
absent for high loading rates and intense loading. 

Qualitative explanation of such a discrepancy can be seen if it is assumed (and it is 
confirmed experimentally) that fracture occurs not at the crack tip itself, but in a certain 
zone ahead of the crack tip because of a complex process of appearance of micro-cracks, 
their mergence, and interaction. The existence of such a zone was clearly demonstrated on 
a brittle acrylic plastic [37], 

Ravi-Chandar and Knauss [38] had observed with high speed micrography that the 
fractured surface formed upon propagation of a crack with a high velocity has three different 
regions: mirror-smooth, matt, and hackle zones. Upon passing the hackle, zone, the crack 
starts branching. The mirror-smooth zone is characterized by a smooth light-reflecting 
surface. In the matt zone the, surface becomes more coarse and in the last hackle zone it 
becomes rough. It is observed that with an increase in loading the size of the mirror-smooth 
zone decreases, whereas those of the matt and hackle zones increase. On the contrary, a 

decrease in the loading level increases the mirror-smooth zone and decreases the matt and 
hackle zones. 

When a single crack propagates in the mirror-smooth zone, its behaviour only slightly 
deviates from the quasi-static growth. In the matt zone, a whole ensemble of cracks 
propagates simultaneously. In the hackle zone, the crack propagation is still of the same 
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physical nature but the zone of micro-cracking becomes substantially larger. Thus, it is 
possible to state that crack propagation under high stresses is controlled by the development 
of micro-cavities and micro-cracks, their accumulation, mergence, and interaction [38]. 

Under conditions of low stress intensity, the crack front exhibits ‘thumbnail’ curvature 
similar to the crack profiles associated with quasi-static fracture processes. At increasing 
stress intensity the crack front straightens out and is identifiable as a front of multiple 
micro-fractures [39]. In the ‘mirror’ zone, crack propagates in the plate with a curved 
front. In the ‘mist’ and ‘hackle’ zones, the ensemble crack front is no longer curved as in 
the ‘mirror’ zone. 

Parton [37] experimentally observed that if crack velocity exceeds a certain value, then 
the arrival of waves reflected from the free boundaries to the crack tip does not affect its 
velocity, implying that the zone of process occurrence has a certain inertia, i.e., it shows 
some resistance to the ‘attempts’ to change its velocity. Based on the boundary element 
analysis in this study, it seems to the author that beyond a certain crack velocity, crack 
behaviour undergoes a distinct change or inertial terms start dominating the square-root 
singular terms. The possibility that the extent of the stress intensity factor during dynamic 
crack growth under transient conditions is more limited than a steady-state analysis would 
indicate was first suggested by the analysis of Ma and Freund [40]. Under severe transient 
conditions, the local crack tip field is not accurately represented by a stress intensity factor 
field. At higher crack speeds, the higher-order terms in the near-field asymptotic expansion 
can have a significant influence on the crack tip stress field [41] even for the case of constant 
crack speed, following fracture initiation. 

During the early stages of crack growth, the other terms may dominate the square 
root singular term under certain conditions. During the early phase of crack growth, the 
transient nature of stress field prevents a complete stress intensity factor from becoming 
established outside the near-tip three-dimensional zone. This feature erodes the value of 
the stress intensity factor concept in characterizing the fracture resistance of the material. 
Though a unique plane stress intensity factor can be identified in both the data and the 
model at each instant of time, even though no fully developed stress intensity factor field 
is apparent. However, this stress intensity factor is not strictly relevant in the spirit of 
fracture toughness testing. 

A long-standing issue of fundamental importance in dynamic fracture research is the 
connection between the dynamic fracture toughness and the crack tip velocity. The debate, 
for the most part, has centered around the question of uniqueness of a relationship between 
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Kf and v. Kobayashi and Dally [42], Rosakis et, al. [43] and Zolmder and Rosakis [44], 
among others, provided data sets that seem to indicate that Kf —■ v relation is ieasonably 
viewed as a material property. On the other hand, the results of Kobayashi and Mall 
[45] and Ravi-Chandar and Knauss [39], based on photo-elasticity and the. method of 
caustics in transmission, respectively, suggest that there is no such correspondence. The 
fundamental difficulty of achieving a well-developed stress intensity factor field may be 
playing a significant role in these differences in results. 

For a crack under suddenly applied uniform crack face pressure, Ma and Freund [40] 
found that for a stationary crack, stress intensity factor field is found only for points closer 
to the crack tip than about 5-10 percent of the distance to the cylindrical shear wave front, 
implying that the extent of the K field is so small at the onset of growth. They found that 
measured and calculated stress intensity factor histories compare very well for relatively low 
crack face pressures, but there is significant disagreement beyond crack growth initiation 
for higher pressures. The time required, before the singular stress becomes a good estimate 
of the total stress, decreases as the delay time increases. Delay time is the time between 
the onset of crack growth and the time instant when pressure is suddenly applied on the 
crack face. The problem may also lie in the assumption of linear elastic mat erial response. 
Both the analysis and the models that underlie the interpretations of those experiments are 
two-dimensional in nature. In either case, the three-dimensional effects that are overlooked 
may have too much of an influence to be neglected in considering the process [40]. 

Conditions of Kf dominance exist either extremely close to the crack-tip or are even- 
tually established at long times after crack initiation. However, the higher order transient 
asymptotic representation is more successful in describing the actual field even at times 
close to the event of crack initiation or at distances relatively far away from the moving 
crack tip. The transient effect is manifested through the time derivative of the dynamic 
stress intensity factor even if the crack-tip speed is constant. 

It is necessary to apply the higher order asymptotic expansion which includes the 
transient history of the crack growth to describe the near tip deformation fields. Also 
the crack growth is indeed controlled by a material related criterion. This criterion gives 
the unique relationship between the dynamic fracture toughness Kf c and the crack-tip 
speed v. The existence of such a criterion was seen validated in a simulation by Liu and 
Rosakis [46] by using the higher order transient expansion, while the lack of the uniqueness 
of a relationship between Kf c and v has been observed when the Kf - dominant assumption 
or the steady state higher order expansion is used [39,45], 
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In the present study, the author has investigated if it is possible to capture this loss of 
square-root dominance by the boundary element analysis. SIF is calculated from crack- tip 
or crack front tractions in BEM. Value of traction at a node on crack front exhibits an 
integrated influence of level of singularity over the elements adjoining it. Any weakening 
of singularity, for example at crack vertex, or loss of dominance of square-root terms at 
high speeds of crack propagation gets reflected in the value of SIF calculated from nodal 
tractions. 

In the past three decades, the application of the finite element method (FEM) to 
dynamic crack propagation has made important advances. Interesting reviews can be 
found in the papers by Aoki et al. [47] and by Atluri and Nishioka [48]. Two different 
approaches followed are based on one of the following two concepts: a stationary mesh, 
which includes a ‘node release’ mechanism, or a moving mesh which normally includes a 
singular element. 

The BEM has appeared as an alternative for elastic fracture mechanics problems. This 
method seems to be a better choice than FEM for elastodynamic fracture mechanics be- 
cause the discretization is restricted to the boundary surface and the concept of singular 
element is simplified. In particular when dealing with dynamic crack propagation, the 
remeshing is conceptually much simpler in BEM than in any domain technique. Classical 
methods like finite difference method or FEM are difficult to apply to the study of crack 
growth because of their intrinsic numerical dissipation of high frequency elastic waves. 

1.3 Advantages of BEM 

The BEM is a numerical technique based on integral equation formulations of the con- 
tinuum mechanics problems. There are two types of integral equation formulations. One 
contains, as basic unknowns, quantities with a clear physical meaning and in terms of 
which boundary conditions are given. This type is called the direct formulation. In the 
second type, known as indirect formulation, the basic unknown quantities have no physical 
meaning. The physical variables are obtained from these sources. In this thesis, numerical 
approach used is based on direct formulation. 

There are some general characteristics of the BEM which represent clear advantages for 
the analysis of static and dynamic continuum mechanics problems. First of all, the prob- 
lem is formulated on the boundary. Therefore, only the boundary has to be discretized as 
opposed to domain techniques, such as the Finite Element Method (FEM) and the Finite 
Difference Method (FDM) which require the discretization of the domain. As a conse- 
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quence, the resulting system of equations is smallei in the BEM. Another con sequence 
is that the mesh generation process, which only involves the surface, is simpler than in 
domain type methods. _ This advantage is important for situations where the geometry 
changes throughout the solution process, for example, crack propagation studies. But the 
final system of equations is non-symmetric and fully-populated. This fact may lead to 
longer computer times, as compared to finite element solutions, in some types of problems. 
A second characteristic of the BEM is that it produces highly accurate solutions on the 
boundary and in particular at any selected internal point. This feature makes the methods 
very appropriate for problems where very high accuracy is required. Among the disadvan- 
tages of the BEM, its difficulty in dealing with non-linear material properties is prominent. 
In such cases the integral equations include domain integrals and the advantages of a 
boundary formulation vanish if extensive zones of non-linear material exist. 

There is an additional important difference between the BEM and tin* domain type nu- 
merical methods. When dealing with infinite or semi-infinite regions, the domain numerical 
methods require a discretization which should extend towards infinity. Obviously, the dis- 
cretization has to be finished at a certain distance where an artificial boundary of some 
kind is located. These boundaries can be used without problems in statics. In dynamics, 
however, spurious reflections of waves which distort the solution of the problem take place 
at the artificial boundaries. On the contrary, the BEM is based on the integral equation 
formulation, which in the case of external regions, consist of integrals extending only over 
the internal boundaries. Therefore, those are the only boundaries to be discretized. The 
behaviour of the boundless domain is ‘naturally’ represented by those integrals over the 
internal boundaries. 

The BEM has special suitability for problems of fracture dynamics. Its capability to 
model steep stress gradients is the main motivation for its use in this study to analyse free 
surface effect for through cracks for both static and dynamic conditions. 

1.4 Present work and layout of the thesis 

The objective of the present work is to investigate the effect of free surfaces on fracture 
characteristics for elastostatic and elastodynamic problems using three-dimensional bound- 
ary element method. Free surface effect is studied for a through center cracked finite body 
for three situations: quasi-static conditions, stationary crack under impact loading and 

running crack under impact loading. The thesis has six chapters and each of them is a 
self-contained unit. 
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In Chapter 2, the three-dimensional direct time-domain displacement formulation is 
described. Software developed on the basis of this formulation has been used to get the 
results of the next chapters. 

In Chapter 3, the effect of free surfaces on fracture characteristics for elastostatic prob- 
lems is discussed. Stress intensity factors are calculated along the crack front for a through 
center crack specimen using nodal traction values at the crack front. Correspondence of 
these stress intensity factors with values of Jf-integral is studied all along the crack front. 
Free surface retardation of crack front is studied for incipient self-similar crack growth and 
critical intersection angle of the crack front with the free surface is estimated. 

In Chapter 4, the effect of free surfaces on fracture characteristics for a stationary crack 
under impact loading is discussed. The variation of dynamic stress intensity factors with 
time is studied along the crack front. The effect of waves travelling in thickness direction is 
probed. Also free surface retardation of crack front is studied for incipent self-similar crack 
growth and critical intersection angle of the crack front with the free surface is estimated 
under step loading. 

In Chapter 5, the effect of free surfaces on fracture characteristics for a running crack 
under impact loading is discussed. Formulation of a moving singular element for three- 
dimensional modelling is described. Variation of dynamic stress intensity factors is studied 
for low and high crack speeds. Free surface effect on DSIF variation at crack vertex and 
other positions on the crack front is studied for straight and curved crack fronts. Also 
critical intersection angles are estimated. 

Chapter 6 presents the main conclusions of this study and guidelines for future work. 



Chapter 2 

TIME DOMAIN FORMULATION 


This chapter describes the direct time-domain displacement formulation used in this the- 
sis. First of all, the general equations governing linear elastodynamics are briefly described. 
Thereafter, boundary integral formulation is given along with detailed numerical imple- 
mentation procedures. Temporal and spatial integration of terms involving displacement 
and traction kernels are described in detail. 

The detailed formulation is given for the sake of completeness. All peculiarities of the 
method used have been brought out in the description. Time domain formulation is a well- 
developed formulation and its details can be seen in several texts [49,50]. In this thesis, 
a new partitioning scheme is proposed for spatial integration of elastodynamic kernels. 
Algorithm for calculating principal value integrals suggested by Guiggiani and Gigante [7] 
is extended to elastodynamic analysis. Also, combination of projection technique with 
semi-analytical integration, used for nearly singular integrals for static problems by Cruse 
and Aithal [8] and Mi and Aliabadi [9], is extended to elastodynamic problems. 

2.1 Equations of three-dimensional elastodynamics 

The basic equations of linear elastodynamics are: 

1. the equations of motion 

&ij,j T — piii (2.1) 


2. the kinematical relationships 
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3. the constitutive law 

Oij - Xfkkkj + 2/u,,- (2.3) 

Also the balance of angular momentum implies that the stress tensor is symmetric = 
c 7 ji ). In the above, Ui(x,t) are the components of displacement, vector at point, x and time 
t, aij and dj are the components of stress and strain tensors respectively, and 6,- is the body 
force per unit mass. Furthermore, A and p, are the Lame’s const, ants and p is the mass 
density. Kronecker’s delta Sy is equal to unity if i = j and zero if i i- j. The Cartesian 
coordinate system (i,j = 1,2,3) is used here. Light-faced letters stand for scalars while 
letters in boldface denote vectors and second-order tensors. 

Above equations can be combined to give the Navier-Cauehy equations 


(A + n) u jtji -1- puijj + phi - friij (2.4) 

which are the governing equations for an elastic, isotropic, and homogeneous solid of volume 
V and surface S. By defining the propagation velocities of the pressure (dilatational) and 
shear (rotational) waves in the solid as c\ = (A + 2 p)/p and cij = p/p, respectively, 


equation(2.4) can be written as 

( C 1 — c 2 ) u j,ij T C‘-2 u i,jj + b{ = Ui (2.5) 

Equation of motion is accompanied by initial conditions 

Ui{x, 0) = 'Uj 0 (2.6) 

Ui(x,0) = v i0 (2.7) 

in V and boundary conditions 

Ui{x,t) = Ui, x e Sr (2.8) 

ti{x,t) = crijTij = I) , x 6 S '2 (2-9) 


on S, where U are the tractions, n^x) is the outward-pointing normal vector and Ui and 
Tj are prescribed values for displacements and tractions respectively. 

It is assumed that Uj and its derivatives up to second order are continuous except at 

propagating wavefronts, where only m is continuous. On these fronts, both the kinematic 
jump condition 


ut\ = -cn k [u iA 


( 2 . 10 ) 
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and the momentum condition 

[U] = - pc[ui } ( 2 . 11 ) 

must be satisfied [51]. In the above two equations, brackets denote difference in the enclosed 
quantity between front and back surfaces of the discontinuous front travelling with velocity 

c. 


2.2 Boundary integral formulations in elastodynamics 

2.2.1 Introduction 

The mathematical basis for all direct integral equation formulations in elastodynamics is 
the dynamic reciprocal theorem. It is not, however, the only method; the same result can be 
accomplished via Green’s function method, variational principles, or weighted residuals. A 
system of constraint equations between the boundary quantities is generated by moving the 
field point about the surface. In order to solve engineering- type problems, discretization 
of the surface is necessary. Subsequent numerical processing of the boundary integral 
equations gives rise to the Boundary Element Method (BEM). 

2.2.2 Dynamic Reciprocal Theorem 

The reciprocal theorem in dynamics specifies a relationship between a pair of elastodynamic 
states. It is essentially the dynamic extension of the classical reciprocal theorem of Beti- 
Rayleigh in elastostatics, augmented by the inertia forces. This theorem was stated by 
Graffi [52] in its present form for bounded regions using Laplace transformations. Rigorous 
proofs for unbounded regions were given by Wheeler and Sternberg [53]. The reciprocal 
theorem first specifies a regular region V in the sense of Kellogg [54] with boundary S and 
material properties p, C\ and Ci- Consider two distinct elastodynamic states A = {v,i,ti,bi\ 
and B = [uj,fj,6j] defined in that region and with initial conditions 
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(2.12) 

^(x,0) = v i0 (x) 


(2.13) 

0) — ^10 (®) 


(2.14) 

'U'i 0) — ^io 


(2.15) 

Then, for t > 0 

1 t{ ^ ZLi&S T / P ^ i H - Viotli + 'U'io'U'io'} ~ =: 
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+ / tj * VfdS 


(2.16) 


The operation * above denotes time convolution, i.e., 


f*g = [ f(x,t — r)g(x,T)dT — f{x,r)g{x,1 ■ r)dr ( 

Jo Jo 


for t > 0 and for two functions / and g. 


2.2.3 Fundamental solutions 


The fundamental solution at point x (receiver) for displacement component i due to a unit 
impulse vector at point y (source) at time r in direction j for a three-dimensional elastic 
solid of infinite extent is solution of equation of motion with a, body force 


bj = 5{x - y)8{T - r) Gj 


( 2 . 18 ) 


where 5 denotes the Dirac delta function and e.j is a constant unit vector in direction j. 
The fundamental solution for displacement at point x and time T is given by 

l r r l/r - 

Gij(x, T;y,r) = — pay - by) J A 8(v - Xr)dX 


+aij {(l/c‘i)5(n - r/ci) - (l/cij)tf(« -- r /<■>)' 
+{bij/cl)5{v - r/c 2 ) 


+(bij/c 2 2 )5{v - r/c 2 )j (2.19) 

where v = T-r; ay — hihj/r 3 ; by = 8y/r\ hi = Xi~ yi and r is the. distance between 
points x and y. The fundamental solution for tractions is given by 

1 r f l / c 2 

H^x, T] y, r) = — ^ - 6cj (5a y - by) / A<5 (u - Ar) dA + (12a,, ~ 2by) { 6 (v - v/c 2 ) 

J l/Oi > 

— ( c 2 /ci) 8 (v — r/ci) j + 2r Oy/ea jd (u - r/c^) — (ra/cj)' 1 A' (u r/cj) 

— Cy (l — 2c^/c?) | A (n — 7'/ ci) + (r/c'i) A (n — 7'/ ^ 

-dyU(v-r/c 2 ) + (r/c 2 ) 8 (v~r/c 2 ))] ( 2.2 


where u T r, hi Xi yy ay — hihjh m n m /r 5 ; cy =-hjiii/r 3 ; dy — (hiUj + 

$ijh m n m )/r 3 ; by = Cy + dy. Also n* are the components of outward pointing normal 
vector and h m n m — h\Tt\ + h 2 n 2 + h^n^. 

The above fundamental solutions obey the causality condition 


Gy- (aj, T, y, r) — 0 if C\(T ~ t) < r 


(2.21) 
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and have the following time translation property: 

Gij (x, T + T 0 ,y,r + T 0 ) = G i3 (x, T, y, r) (2.22) 

Love [55] concluded that Stoke’s solution yields correct results only when the input quan- 
tities are continuous at the propagating wave front. In this case, both kinematic jump 
condition and momentum conditions are satisfied as well. This implies that if the input 
excitation is a step loading, it must be modelled as a ramp loading during the first time 
step. Also very small time steps should not be used because they result in non-vanishing 
dilatations and rotations at the wave front. 


2.3 Boundary Integral representation 


By utilising the dynamic reciprocal theorem with the untilded state being the real one and 
the tilded state being the fundamental singular solution, Love’s [55] integral identity may 
be written in the form 


Cij(y)ui(y,T ) = J [ Gij*U(x,T ) - H i3 * Ui(x,T)]dS(x) + p j G i3 * bi(x,T)dV(x) 

GijV i 0 (x ) + GijU io (x ) dV(x) (2.23) 


+P 


In the above, x and y are receiver and source points respectively, and u i0 and Vi 0 are the 
initial displacement and velocity, respectively. The free-term tensor components Cij can 
be expressed as 

Cij — 5ij — 7 ij (2.24) 

where is the Kronecker’s delta and 7 ^ is a discontinuity or jump term. 

C-matrix depends on the local characteristic surface of the boundary point with respect 
to the domain. It depends on the solid angle of the boundary point. Mantic [56] has given 
a general closed analytical formula of the C- matrix for the case of any finite number 
of tangent planes to the boundary of the body at a non-smooth boundary point. For a 
smooth surface, value of Cij is .5 An indirect method is also often used to calculate 
C-matrix as illustrated by Cruse [57]. C-matrix remains the same for both elastostatic 
and elastodynamic problems because it depends on the local geometry and Poisson’s ratio 
of the material. 

The above integral representations are exact statements of the problem that was origi- 
nally cast in the form of partial differential equations. No approximation has been made, 
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and solutions of the integral representations satisfy the original differential equations as 
well. 

In this work, elastodynamic problems, having zero initial conditions and zero body 
forces, have been considered. For zero initial conditions and zero body forces, the boundary 
integral formulation for transient elastodynamics reduces t o 

Cij (; y)u { (: y,T ) = J [G ij (x,y,T)*t i {x,T) - Il i:i {x,y,T) * ■u,(xJ')}dS{x) (2.25) 

where 


Gij * ti 


Hij * Uj 


rT 


Gij{x,T]y,T)ti (as,r)dr 
Hij (*, T ; y, r) u t {x, r) <1 ; 


(2.26) 

(2.27) 


2.4 Numerical treatment of boundary integral equa- 
tions 

For zero initial conditions and zero body forces, the integral statement. (2.25) expresses 
the unknown variable in terms of boundary variables. The solution of such an integral 
equation is not possible, in general, except if it is done numerically. The geometry of the 
problem and the field variables (displacements and tractions) must, be discretized to solve 
non-trivial problems. 

The time-domain BEM in three-dimensions consists of two basic steps as explained 
in [49,58,59]: 

1. A discretization of the real time axis into a sequence of equally spaced time intervals 
with the assumption of linear variation of displacement, and constant, variation of 
tractions over each time interval. When piecewise linear variation is used for both 
displacements and tractions, the time-stepping scheme becomes unstable [GO, Cl] 

2. A discretization of boundary into 8-noded elements. 

On the basis of these discretizations a time stepping solution of equation (2.25) can be 

established for the boundary displacements and tractions over each element, and for each 
time step. 
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2.4.1 Spatial discretization 

The boundary is divided into a number of elements. The geometry, in the present work, 
has been described’ by continuous quadratic, eight-noded elements. Triangular elements 
are treated as collapsed quadrilateral elements. The cartesian coordinates of a point on 
the surface of a boundary element are given in terms of the element nodal coordinate X ia 
as 

Xi = M a (j]k) X ia (2.28) 

where i = 1,2,3; k = 1,2; and a = 1, 2, . . . , 8. 

In the above, M a are the shape functions defined in the intrinsic or local coordinate 
system %. The shape functions for the quadrilateral element are: 

f -25(1 + fo)(l +%)(& + 77o-l) if o=l, 3, 5, 7 
M a = < .5(1 + £ 2 )(1 — tjo) if 0 = 2,6 (2.29) 

[ .5(1 + £ 0 )(1 “ V 2 ) if o = 4,8 

where £ 0 = £.£ Q and rj Q = 7].T] a with £ and rj being the two independent coordinates and 

(£q, Va) the coordinates of node a. 

The transformation from the cartesian coordinate system to an element intrinsic coor- 
dinate system is through the Jacobian matrix given by 

Jij = ^X ia (2.30) 

dVj 

2.4.2 Temporal and spatial interpolation of displacement and 
traction 

The time span of interest is divided into N equal time steps of size AT so that 

T n = nAT, n = 1,2,..., AT. (2.31) 

The boundary variables are interpolated in terms of their values at nodes by 

N 

u,(x,T) = [Ml <■'(*) + «?(*)] (2.32) 

71=1 

N 

ti(x,T) = 2 bf tr l (x) + & tfW] (2.33) 

n — 1 
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and 


u"(x) = £> r K(*K r . 

C=1 

r*(«) = Ewifwi”!""". 


1 . 3.3 


/■ 1,1?. R 


(2.34) 

(2.35) 


C=1 


where u™ and u\ n represent nodal values ol displacement at time nodes n and (/?. — 1) 
respectively. Similar interpolation is done for traction variables also, lb-re. .V is the total 
number of time steps. Mj and M F are the temporal interpolation functions, related to 
local time nodes I (initial) and F (final), and are of the following, form for pica-wise linear 
variation: 


where 


<£n( T ) = 


or 


II 

vT 

T _ r 

V Mr] 

( 2.3(1) 

Mp = 

1 «•> 

(2.37) 

f 1 if 

(n - 1)A7’ < r < n.\'T 


\o 

otherwise 


H (t- 

(n-l)A'J’) - II [r u.\D 

(2.38) 


H being the Heaviside function. 

The displacements are interpolated by piecewise linear functions. Tractions an- repre- 
sented by piecewise constant functions which allow for sudden jump of their values from 
step to step. It has been shown, in the case of BEM bast'd on displacement integral 
equation, that this is the best interpolation scheme to linear order. {f>< L E I i For tractions, 
temporal interpolation functions are of the form 


M/ = 0.5 cl> n (r) (2.39) 

Mf = 0.5 (j> n {r) (2,10) 

For spatial interpolation of boundary variables, an approach |(i2| which combines the 
advantages of continuous and discontinuous elements in a single formulation is used in 
this paper. In this formulation, expressions of the shape and interpolation functions are 
unified in a form suitable for fully continuous, fully discontinuous or transition elements 
(i.e., discontinuous only along one or more edges). Three types of boundary nodes are 
identified depending on the function they fulfil: 
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1. Geometrical nodes: those used for interpolating the boundary geometry. 

2. Interpolation nodes: the points on the boundary which define the interpolation ( 
the variables.' The system nodal values correspond to these nodes. 

3. Collocation nodes: where the boundary integral equation is enforced. 

These three kinds of nodes coincide at the same point in the commonest formulation c 
the boundary elements. One property of the boundary element formulations is that th 
discontinuities of the variable across the element edges do not invalidate the convergenc 
of the technique [63]. Collocation and interpolation nodes for continuous elements coincid 
with the geometrical nodes. Collocation nodes and interpolation nodes for discontinuou 
elements and discontinuous side of transition elements are situated inside the elemeni 
For quarter-point element, geometrical and interpolation nodes coincide, but collocatio 
corresponding to nodes on crack front is done at nodes situated inside the elements o 
either side of the crack front. To reduce the number of equations and enhance accuracy 
equations corresponding to collocation nodes before and after a given interpolation nod 
are added up [64], 

Interpolation functions 

To determine these interpolation functions, it is necessary to define the position of the nods 
points within the element, which can be done by defining a certain number of parameter 
based on the local coordinate system of the element. These parameters for the quadrilaters 
elements are as follows: Six parameters are needed to define the element nodal positions 
ax, a 2 , 03 , a 4, 05 and a 6 , with a, = 1 — 6 i; [i = 1 , 2 , 3, 4), a 5 = 02 + a 4 , and 06 = Ox + 03 (Fig 
2.1). The values = 6 2 = 63 = 64 = 1/3 are assumed for quadratic interpolation. Th< 


interpolation functions for eight-noded elements are as follows: 

M{ 2) = (o 2 - 0 (a 4 + 0 (03 - v) / (02 o 4 o 6 ) (2.41 

M( 4) = (a 3 - 77 ) (ax + 77 ) (o 4 + £) / (ax a 3 a5) (2.42 

M( 6 ) = (o 2 - 0 (o 4 + f) (ax + 77 ) / (a 2 a 4 a 6 ) (2.43 

M( 8 ) = (a 3 - 77 ) (ax + 77 ) (a 2 - f) / (ax a 3 a5). (2.44 

M(l) = (a 2 - £) (a 3 - 77 ) /a 5 /a 6 - M(2)a 2 /a 5 - M( 8 )a 3 /a 6 (2.45 

M( 3) =• ( 03 - 77 ) (a 4 + £)/a 5 /a 6 - a 3 /a 6 M (4) - a 4 /a 5 M ( 2 ) (2.46 

M( 5) = (01 + 77 ) (04 + 0 /^ 5/06 - Ox/o6 M (4) - a 4 /o 5 M(6) (2.47 

M( 7) = (a 2 - 0 (ax + 77 ) /o 5 /a 6 - ax/a 6 M ( 8 ) - a 2 /o 5 M ( 6 ) (2.48 
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7 6 5 



Figure 2.1: Q uad ril at oral oiomonl 


2.4.3 Time- marching scheme 

The displacement boundary integral equation can be written as 
rT " r i-r N . , r 


HijUi] dSdr 


[(’ijl'i //,/»,] rfSdr ( 2 . 49 ) 


CijUi(y,Tj\r) — / I [Gijti 

J Tn-Js 

where the integral on the right hand side is the contribution duo to past dynamic 

This equation is also an exact statement since no approximation has yet been introduced. 

However, in order to solve the above equation, one has to approximate the variation 

of the field quantities, in addition to the usual approximation of spatial variation. The 

time integration in the above equation is done analytically, and the spatial integration is 

one numerically. A system of algebraic equations is obtained at. time step /V through 

nodal collocation, l.e„ by allowing y to coincide sequentially with all the nodal points of 
tne boundary. 

The boundary integral equation for the first time step is 

Cij u i(,V, Ty ) ~ ^ [ Gijk - HijUi] dS dr - () 

After integrations, the resulting system of equations lias the fori. 


( 2 . 50 ) 


rm 


{t 1 } - [J3j] {„■} + [G>] {(«} - [Hj] {„»} = {0) 


( 2 . 51 ) 
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where [G] and [H] are coefficient matrices resulting from the surface integrations, {f 
and {«} are the tractions and displacements at all nodal points, and superscripts denot 
the time steps. For [G] and [ H ] the superscript denotes the time step at which they ar 
calculated, and the subscript denotes the local time node (J or F ) during that time interval 
After the usual assembly process, the resulting system of equations has the form 

PM {V 1 } - [sy {y 1 } + K] {*»} - [b‘] {r»} = {0} (2.52 

where [A] and [B] are the matrices related to the unknown and known field quantitie 
respectively. {X} and {T} are vectors of unknown and known field quantities respectively 
Since all the unknowns at time T — 0 are assumed to be zero, equation reduces to 

PM {A 1 } = [By {V 1 } + \B}\ {y°} (2.53 

Now, the boundary integral equation for the second time step is considered. 

CijUi{y,T 2 )- f [ [GijU - HijUildSdr = f [ [Gijti - Hiju^dSdr (2.54 
Jti Js Jo Js 

If the time difference (T 2 — Ti) is the same as (7\ — To), the resulting coefficient matrice; 
of the left hand sides of equations (2.50) and (2.54) become identical. This is due to tim< 
translation property of fundamental solutions, G i; - and Hij, which contain time function; 
with arguments (T — r) and therefore, the convoluted integral corresponding to the interva 
Ti < t < with T = T 2 is identical to that of the interval T 0 < r < T\ with T = T\. 

The right hand side of equation (2.54) is evaluated at time T = T 2 with the timi 
integration over the interval To to Ti and thus, provides the effect of the dynamic histon 
of the first time interval on the current time node, i.e., T 2 . 

Now, the resulting system of equations for the time node T 2 is of the form 

[A' F }{X 2 }-{B' F ]{Y>}+[G}}{t'}-[H}]{u'} = -[& F ]{t'}+[H* F ]{u'} 

- [Gj] {t 0 } + [. H ]} {u 0 } 2.55; 

Equation can be rearranged such that 

PM {A' 2 } = [By {y 2 } - [g; + cy { t q + p; + b 2 ] {«•} + [«?] {«»} - [g?] py 

'2.56 

In the above equation, all the quantities on the right hand side are known. Therefore, th< 
unknown vector {A!} at time T 2 is obtained by solving the above system of equations. 
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Similarly, for the N th time step, the boundary integral equation is written in a dis- 
cretized form as given below: 

[4] {*»} - [j*] { r*} = [«f 1 "" 'll" 

+ tH7lK}-l«Vl{'“} (2- 57 ) 

or 

[4]{^} = [b;..]{v"} + {r k } cm) 

where vector { R N } is the effect of the past dynamic history on (.lie current time node. The 
above equations can be solved to find the unknown vector {A' ,v } at time Y'.y. 


2.4.4 Analytical temporal integration 

For linear time interpolation, field variables are expressed as 

N 

/(*, r) = £ K7 n ” 1 (*) + Mjr M | (2.59) 

71=1 

where f n {x) represents the spatial variation of the field variable /(:r, r) at i ime T n ( nAT ), 
and Mj and Mp are the temporal shape functions defined earlier. 

The transient kernels, involved in boundary integral equations for clast odynamie prob- 
lems, have one or more of the following time functions embedded in them: 


Time function (1) : 5 (T - r ~ r/c) 

Time function (2) : J^ 2 A 6{T -r ~ r/c ) dA 

Time function (3) : 5 (T - r - r/c) 

Time function (4) : S (T - r - r/c ) 

Time function (5) : 5 (T - r ~ r/c) 

where c is either the pressure wave velocity Cj or the shear wave velocity »•->, and rf is the 
Dirac delta function. 

Using equation (2.59), the time integrals related to any of the above listed functions 
can be expressed as 


f g{T — t — r/c) f (x, t) dr 
J o 



E r-\*) [ 

n=l L «'(» 


nAT 

(n-l)AT 


1 r/c) f (x, r ) dr 
M } 1 (r) g (T - r - 


r/c ) dr 


i 
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nAT 


+ f n {x) Mp (t) g (T — t — r / c) dr 

' (n-l)AT 


(2.60; 


An important characteristic of the transient dynamic kernels is the time translation prop 
erty. Because of this property, at each time step, only the effect of the dynamic history 
of the first time step on the current time node needs to be evaluated, i.e. , at each tim< 
step, the analytical time integration has to be done only for n = 1. Thus, equation (2.60 
reduces to 

f Ar r AT 

J o g{T -t - r/c) f (x,r)dr = /° (*) y (i - ~£rpj 4>i i T ) 9 {T - t - r /c) f (a?,r) 

+ Pin) J (at) ^1 ( r )‘ 9 ( T ~t _r / C ) f (*' T )^' e - ( 


The time integrals on the right hand side of above equation can be rearranged into th< 
following two types of integrals: 


• Type A: f AT fa (r) g (T - r - r/c ) dr 

• Type B: f AT rfa (r) g (T - r - r/c ) dr 


The convoluted time integrals of type ‘A’ and type ‘B’ with the time functions (1) and (2 
are evaluated analytically as follows: 


Time function 1 
Type A 

AT 

fa (t) 5 (T - r - r/c ) dr = fa (T - r /c) = H {T - r /c) - H {T - AT - r/c ) (2.62 
Type B 

AT 

rfa ( t)6{T-t - r/c ) dr = (T - r/c ) fa ( T - r/c) (2.63; 




Time function 2 
Type A 


AT 



A 5 (T — t — X r) fa (r) dAdr 



A (f> (T - Ar) dA 
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<Ai ('/' \> ) 


1 /<■;• 


II (T Ar) 


i /.-i 


A' 2 


— II {T A 7’- Ar) 


1 /<’ 


( 2 . 64 ) 


where 

1/C2 
l/ci 

The second term on the right hand side of equation (2.64) can bo obtained in a similar 
manner by replacing T by T — AT in equation (2.65). 




r 0 if T < r/r, 

S - Sf » ’'A'. < T < ’■/<'•-' (2.65) 

1 sj - s? if r> 


Type B 


AT /* 1 /C2 



0 Jl/Cl 


A 5 (T — r — Ar) r(j>i (r) dAdr 


f 1 ‘ \(T~ \r)<f) X (7* Ar) <1A 
l/ci 

1./C2 r 1 / <’2 

TA0, ('I 1 - Ar) (1A • / rA'V’i (7’ - Ar) dA 

1/ci !/n 

l/f-2 

( 2 . 66 ) 



[A 2 1 

i In 

r A :i 

- 

T 

J<h ( 1 ' - Ar) 

- r 

i/n 

a 

• At) 


i/«i 


The terms on the right hand side of equation (2.66) are evaluated in a similar manner as 
that of equation (2.64). 


Time function (3) 


enAT 


'(n-l)AT 


S(T-r- r/c) f(x,r) dr = (T .. 


( 2 . 67 ) 


Time function (4) 

The temporal integration involving the time function (4), i.e., S (T ~~ t - r/c), is approxi- 
mated by using a backward finite difference scheme as follows: 


r»nAT 


/(n-l)AT HT ~ r - T ^ f (*’ 0 iT = ~ (X) + ~ (T ~ T/C) (2 ' 68) 
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Time function (5) 
Similarly, 


r T 5 (T — r - r/c) / (», r) dr - M (T - , 

(n-l)AT (AT) 

(2.69 


If temporal integration results are closely examined, a particular radial variation is seen 
All temporal terms can be expressed in the form a + br + £ after analytical integratioi 
with time as shown below: For the first time interval, 


MJ(t) = 

(l Ar )*« 

(2.70 

ACJ.(t) = 

AT ^ ^ 

(2.71 


Hence, for example, 

r AT 


1 0 


5 (T — t — r/c) M) (t) dr = ^ 1 


_ / (1 - sO + Sf ' r if C ( T - A T)<r<cT 

0 otherwise 


(2.72; 


And 


cAT 


6(T-r- r/c) Mp (r) dr = <j ^^r) r lf C ( T AT ) < r < cT 


otherwise 


(2.73; 


Similarly, 


n l/C2 

Mj (r) 5(t- Xr) AdA = A - 

Id 


B 


AT 


(2.74; 


where 


A = iA 2 7/(T- Ar) - ^A 2 H(T- AT- Ar) 


= £>-£ 


(2.75) 


And 


D 


{ 0 if r > ci T 

i(4-4) if r<C2T 

Hi + 14 « ^<r 


(2.76) 


£ 


< ciT. 

0 if r > ci (T - AT) 

i(i-i) if r<c 2 (T-AT) 

- 1 2 ^+ 1 2 {J ^ if c 2 (T — AT) < r < Cl (T — AT) 


(2.77) 
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Similarly, 

B 

where 


F-G 


T 

( 1 

1 A 

1 1 

1 1 

1 \ 

2 

U“ 

■sJ 

3 1 

U' 



G = 


1 T i R 1 , T 3 

2^ + I ^ + 6^ 


“ nr 


if r > c.\T 


if C'/T < r < c{r 


if r > c, (T - AT) 
if r < c -2 (7’ - A7 1 ) 


(2.81) 

(2.82) 


| Jr + f + ^~&p- ( T + 2AT) if <?2 (7 1 - AT) < r < rq (7’ — AT) 


Similarly, 


/•at 12 

Jo l/<, M ’" {T) XS{1 ~ XV) d/A " "AT (2.84) 

where R is defined in equation(2.79). 

Thus, all temporal terms can be expressed in the form a + hr + A after analytical 
integration with time. Also it is noticed that for the distance reached by both shear and 
compression waves travelling from a source node up to a time AT, coefficient, of 4, is zero. 

2.5 Spatial integration 

In this section, details are given about a new partitioning scheme that is proposed in this 
thesis for spatial integration of elastodynamic kernels. A typical surface element is shown 
m figure 2.2. At first the element is quiescent since no signals from the source point, y have 
yet arrived. Spatial integration of the displacement and traction kernels takes place during 
the time interval t m denoting the arrival of the pressure signal. The area, over which the 
integration takes place is the band between radii c x t m and c 2 (t m - 1). After passage of 
the shear wave, the element is quiet again and there is nothing to integrate in view of the 
causality condition. The area covered by a wave in a time step, starting from a particular 
source node, does not contribute to the integrals in the later time steps. 

The^kernels m elastodynamic case have discontinuous derivatives at points such that 

r/Cl ! 2) = kAt ^ mte S er )- Therefore, the standard numerical methods, like Gaussian inte- 
gration, give poor results if applied over the whole element, say, E { . Instead one has to 
partition each element E, into sub-elements E$ : $ = u with E$ = {(* - l)At < r/c < 

t}. Then, an integral over E { is the sum of sub-integrals over E\, each sub-integral being 
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evaluated using standard numerical methods. The partitioning is reasonably tractable for 
2-D problems. It is clearly a purely geometrical problem (find the intersection points of a 
given curve C with circles centered at a given point y and of equally spaced radii) . On the 
contrary, its three-dimensional equivalent (find the intersection curves of a given surface 
S with spheres centered at a given point y and of equally spaced radii) is extremely com- 
plicated. However, Karabalis [59] had performed the spatial integrations by partitioning 
the receiver elements into very small rectangular sub-elements. With this, the possibil- 
ity of only certain portions of the source element being ‘active’ also could be taken into 
consideration. 

In this thesis, a new scheme is proposed for partitioning of surface elements for 3-D 
problems as follows: 


2.5.1 Partitioning scheme 

Consider a boundary element S p . It is mapped onto a square region R ( Fig. 4.1) in the 
parameter plane having local intrinsic coordinates £ = (£i,£ 2 ). Suppose source point is 
y and a typical receiver point is x. Hence, a typical integral over the surface S p can be 




Figure 2.3: Image in the parameter plane and polar coordinates 

written as 

1 = / f h(x(0,T-,y,T)N‘(i,x)M(T)dT}dSlx) 

JSpJ 0 L 

= [ \f ^(x(e),T;y,r)M(r)drl N^^,x)dS(x) 

Js p U o 

= f [ [^(^(e),r;tj ) r)M(r) f ir]^'(0.;(()<16<16 (2-85) 

J RJ 0 

where Vij is displacement or traction kernel, M(r) is temporal interpolation function, and 
N c (£,x) are spatial interpolation functions of field variables. 

Coordinates of a receiver point on a boundary element are given by a parametric rep- 
resentation in terms of shape function M c (^ 1 ,^ 2 ) and nodal coordinates x'\ 

= * = 1,2,3 (2.86) 

C 

Let us indicate by rj = (%,%), the image of the point z on the element. Point 2 is closest 
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to the source point. The first, second and third order derivatives are obtained by 

dxi dM c „ 


d& d£ k Xl 

d 2 Xi _ d 2 M c 

W k ~~d$ 

d 3 Xi d 3 M c 


X; 


Xi 


(2.81 

(2.85 

(2.8S 


By employing a Taylor expansion about the closest point, it is easy to establish formula 
of the form 


Xi 


Zi + 


+ 


+ 


d xj 

QZi 

d 2 x 


e =» 7 


(6 - m) 


£=77 


3£i 


(6 ~ Pi) 2 

Ok ' 


£=77 


5 (f i ■ ,?l)2(& “ m) 


(6 “ Pi )(6 -%) + 




£=77 


0$ 


(6 - pa) 1 




?=7? 


1 / . \//- \2 ® 
+ 5 «1 -Tl)(&- %) 


?=7) 


(2.90 


Since geometry is interpolated with help of quadratic functions, above formula is exact 
The terms 4^- and do not appear in the above expression as these are equal to zer 
due to quadratic interpolation of geometry. 

Following a common practice in the BEM, polar coordinates (p, 9) centered at rj ar 
defined in the parameter plane as 


6 = V\ + P cos 9 
£2 = V 2 + P sin 9 


(2.91 

(2.92 


so that d £ 2 = p dp d9. Hence, 

^ - Vi = e, + p 


dxi „ dxi 


1 cos0 + — sin# 

L<96 d£ 2 


+P 


d 2 Xi cos 2 0 9 2 Xi 


3 2 :Ci sin 2 9 1 

W~2~ 


2 

1 


+ 




cos 9 sin 9 + 


+p 6 ^ cos 2 9 sin 9- ^ - ^ - - + ^ cos 9 sin 2 9 w Jj% 

[2 d({d& 2 

ej + pAi ( 9 ) + p 2 Bi (9) + p 3 Ci (9) 


d£i <9£l 


(2.93 
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where e* = zi~ yi- Hence, 
r 2 = Zixi-yt ) 2 

= Ee- + 2pE(A t (0)ei) + p 2 (SAf (6>) + 2E(e i J3, ; (6'))) 4- 2 p'\V* A ,{0)11^0) I >-(«•, -(^(0))) 
+p 4 (2S(4(0)C,(0)) + EJ3?(0)) + 2p 5 'E(Bi(9)Ci{0)) + /EC’f(^) (2.94) 

For integration, the integration square (-1 < & < 1) is subdivided into triangles with 
common vertex at r? = (771,772)- The integral (2.85) of a typical kernel over each triangle is 
performed as 

r0'2 ppmax 

/= / / F«(ft0)d(*W (2.M) 

«/ 0 

where 

(p, 0) = Wij (y, x (£ (P, 0))) iV c (* (p, 0)) MCI/'. 0)) (2-96) 

where = J Q r Vy- (as,T; y,r) M(r)dr. The limits 0) and 0 a are angles subtended by the 
arms of the triangle on the vertex rj — (771,%)- Here p max is the value of p on a particular 
ray making angle 9, where it meets the opposite side of the triangle. 

First of all, Gauss integration is performed over 9 direction. Corresponding to the value 
of 9 computed as per Gauss points, p max is calculated for the particular ray. The values 
of p on a particular ray corresponding to the distance of c.^iAt or (\>nAt are computed 
by a combination of bisection and Newton- Raphson methods applied on equation (2.94). 
For each ray, the value of p are found corresponding to positions of heads and tails of the 
wavefronts of compression and shear waves, starting from the source node, for the element 
and time interval under consideration. The ray is divided into three different /.ones for 
separate Gaussian integration. The three different /ones are: 

1. both compression and shear wave cover part of the ray. 

2. only compression wave covers part of the ray. 

3. only shear wave travels that portion of the ray in that particular time interval. 

In a time interval, only those areas contribute to the integral which were not affected 
by the particular waves in previous time intervals. Thus, jump in the kernel at points such 
that r = Ci( 2) nAt is handled accurately. 
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Non-singular and singular integrals 

Analytical integration of the matrix coefficients in equations is not, in general, possible anc 
therefore, numerical quadrature must be used. Thus, fundamental solution-shape functioi 
products are approximated by application of the Gauss- Legendre quadrature formula. [65 

r i k 

/ (2.97 

k=l 

in which r) k and w k are Gaussian points and weights, respectively. Application of th< 
above one-dimensional formula along two orthogonal directions allows integration ova 
quadrilaterals. 

One can distinguish two basic cases of integrand behaviour, non-singular and singular 
In the former case, the distance r between source-point y and receiver (or integration 
point x on the surface of the element is never zero. In the latter case, y coincides with on< 
of the receiver points and r becomes zero. 

2.6 Evaluation of Principal- Value Integral 

Principal value singularity occurs for traction kernel Hij when source element becomes alsc 
a receiver. All temporal terms can be expressed in the form a + br + ^ after analytica 
integration in time, as explained in section 2.4.4 on temporal integration. Also a, b, d are 
constant in the respective zones as explained in section 2.5.1 allowing use of techniques 
developed for integrating elastostatic kernels by splitting the radial integration over three 
separate zones. 

It is well known that integral involving the traction kernel is strongly singular only ir 
the first time step. It is noticed that for the distance reached by both shear and compressior 
waves travelling from a source node up to a time T, coefficient of ^ is zero. Thus, the 
temporal terms are not contributing to the singularity in r. In the later time steps, the 
integral is regular as the region, close to the source node and traversed by both the waves 
up to the previous time-step, does not contribute to the integral. The Cauchy principal 
value integral is treated in the same way, as done by Guiggiani and Gigante [7] for the 
elastostatic case. The non-linear transformation suggested by Doblare and Gracia [66] is 
also incorporated. | 

Elastodynamic displacement boundary integral equation for a three-dimensional do- 
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main Q, bounded by the regular surface [54], can be expressed as 


lim f [Gij(x, T; y, r) * U{x, T ) - #*/(*, 7’; y, r) * u,{x y 7’)] <\S(x) - 0 (2.98) 

+i (s _ e£)+ s £ 

If r = | x _ y| represents the distance between the source point y and receiver point x, 
the kernel Gij shows a weak singularity when x — > y, while the kernel //,, shows a strong 
singularity of the order of 0(l/r 2 ) in three-dimensional problems. 

Since in equation (2.98) both the source point y and the integration point, x lie on 
the surface S, a limiting process is necessary. Since equation (2.98) stems from Green’s 
second identity, it may only be formulated on a domain not including t he singular point 
y. The situation is shown in Fig. 2.4, where a (vanishing) neighbourhood v, of y has been 



Figure 2.4: General vanishing neighbourhood around the source point, 


removed from the original domain Cl. The integration is thus performed on the boundary 
S e — (S - e e ) + s e of the new domain Cl e = Q - v c . The integration must be done before 
taking the limit. Since equation already states that the value of overall limit is zero, all 
divergent parts get cancelled [7,67]. 

It is also not necessary to restrict the shape of v e . It may have any shape, provided y 
is an exterior point for 0 £ , and S ( is a regular surface in the sense of Kellogg [54], 
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However, once the shape of the vanishing neighbourhood v e has been selected, both < 
and s £ are uniquely determined, and their shapes must be preserved while t — > 0 + , i.e 
the maximum chord of v € approaches zero. More importantly, the shape of e £ must t 
consistent with the shape of s £ throughout the process. 

The most convenient shape for s £ , that is a sphere centered at y and of radius e, : 
selected. The selected shape of s e also enforces the shape of e £ , which becomes a symmetri 
neighbourhood on S around the singular point y. Although the value of limit taken as 
whole in equation is completely independent of the selected shape of v e , the value of eac 
term in equation(2.98), taken separately does actually depend on the shape of either e £ c 
s £ . 

Since s e is a sphere, the limits of all integrals on s £ in equation(2.98) can be evaluate 
explicitly. Letting e — > 0, the following boundary integral equation arises 

Cij (y) Ui (y, T) = [ [Gij {x,T-,y,r) *U(x,T) - Hij (x,T]y,r) * Ui(x,T)]dS(x) 

JS-e e 

(2.99 

where the so-called free-term coefficients are given by the limit of the traction integrc 
on the surface s £ of v e , while displacement kernel gives a zero value over s £ . Due to th 
chosen shape of v e (a sphere for three-dimensional modelling) which enforces the symmetr 
with respect to y of the vanishing neighbourhood e £ on S, the surface integral on left- ham 
side of equation(2.99) can be defined in the sense of a Cauchy principal value. 

The set of all adjacent elements connected to the singular collocation point y is definei 
as S p . Those elements belonging to S p are indicated by S™(S P = U m S ™). The boundar 
element S™ is mapped onto a region R m of standard shape in the parametric plane(i.e., ; 
square). 

Now, the geometry of an element is described in terms of shape functions and noda 
coordinates: 

**(?) = * = 1.2,3 (2.100 

c 

Displacement and traction are described by shape functions N c (( 1 ,^ 2 ) of local intrinsic 
coordinates £ = (^ 1 ,^ 2 ), so that 

N 

Ui (x , t) = *^2 + Mpu'lix)} (2.10L 

71=1 

and 

<(*) = * = 1.2.3 (2.102: 
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[«*)]«!" ’ V - ' '■'’•3 (2.103) 

where u™ and uf~ l)c represent nodal values of displacement al time node n and ( n _ n 
respectively. 

Similar interpolation is done for traction variables also. Aft or disrret izat ion of geometry 
and interpolation of boundary variables, discretized version of equation(2.dd) is obtained 


«f'”(T)[iV“(e(!/))Ci J (y) + J s l m*,T;y,r) W?(rMr,W(»)) 

f ({(«)) Af-( rWr ,|.s-(*) 

J Sp JO 

Wij (»,7 , ;|/(0, r) .V' (0 M’/{; hi; | ./ #, (0<I< 


-EE « 

b c (b,c)y£(p,a) 

-EE 



b c (W#(p,a) jRb J ° 



R b J 0 
r 


[%(*//’; ?/(Oo) hi, 

+ ? ? ^ ,n ' 1 L„ l [(?ii (x ’ T; ’ ?/( ^ r) A "'^ )A/ n ? ) (! ; 1 Ai)<\i 

+ Ljo ^(*» T ^(0, r )^ r (0^r)<lT|./ 6 {O«li (2.104) 

::rr \ r “ “"*• ** *«*«• 

.ppmg between the boundary points and the intrinsic coordinates (. the SVIn „ lrt ry of 

the vanishing neighbourhood hidden in any CPV rtetbh.i,,, . , . 

Tn Am, of fr, 1A A „ y v definition must be maintained. 

singular (coBoca ! ’ ’ ““ "* <* a " •*»»«.» connected to the 

o \ r y cement h b m the parameter space £ (i e I lily f_i ll 

square). In the next section, attention is focussed on CPV inLgra, on S)', ' " ' 

Strongly Singular integral on surface element 

(s ’ = aro “ by * 


S„ is 


Is, l N Wf) ’ T ’ V’ T > N ‘ («*)) Af(r)dr] d.9(x) 
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lim 

e ->0 


Sp — 6c 


Hij ( x,T;y,r ) M(r)dr 


N a ({(x))dS(x) 


lin *E f [ T H ij (x,T ] y,r)M(r)dr} N a (x)) dS (x) 

e ~>° ~r Js™-e™ Jo 

m p € J 

lim 2 [ f WuM(),T-,y,T)M(T)<lr]N‘‘(()r'(t)<l£ 

£_>0 m JO 


(2.105) 


where e™ represents the part of e e lying on S™. The plane region R m (of standard shape) 
is the image of the surface element S ™ in the transformed plane (£i , £ 2 ) and accordingly, 
o™ is the image of e™ (Figure (2.5)). It can be noticed, in general, that a ™ does not have 
a circular shape (although its polar radius still tends to zero with e). In equation (2.105) 

, all elements S™ around the pole y must be taken together, as any single contribution is 
unbounded in the limit. Now, elastodynamic kernels contain time varying functions. All 
temporal terms can be expressed in the form a + 6r + ^ after analytical integration with 
time as already explained. Also a, 6, d are constant in the respective zones as explained in 
section 2.5.1, allowing use of techniques developed for integrating elastostatic kernels by 
splitting the radial integration over three separate zones. 

It is well known that integral involving the traction kernel is strongly singular only in 
the first time step. It is noticed that for the distance reached by both shear and compression 
waves travelling from a source node up to a time AT, coefficient of ^ is zero. In the later 
time steps, the integral is regular as the region, close to the source node and traversed by 
both the waves up to the previous time-step, does not contribute to the integral. 

Thus, Cauchy principal value integral for elastodynamic case can be treated in the 
same way, as done by Guiggiani et al [7] for the elastostatic case. Next section shows 
the transformation of the original CPV integral into an element-by-element sum of regular 
integrals, each one expressed in terms of intrinsic coordinates. 


Principal-value calculation 

A representative term of strong singularity present in the traction kernel can be represented 
as ^ 2 . This section shows how the singularity is subtracted out and integrand is regularised. 

A polar coordinate system (p, 9) centered at 77 , image of the singular point y , is intro- 
duced in each R m . 

Define 


6 = + P cos 9 

6 = V 2 + P sin 9 


(2.106) 

(2.107) 




Figure 2.5: Image in the parameter plane of the boundary element and of the vanishing neigh- 
bourhood 


thus, having 


= pdpdO 


( 2 . 108 ) 


The integration on element S™ is now considered separately, provided the limit process 
has not been performed yet (its value would be unbounded). 


I R m -a “ 


n P m {0) 

F™ (p, 0) d pdO 

4 


( 2 . 109 ) 


where 


Fij (P) 0 ) — Hij {y ) & (£ (p) $))) N a (£ (p, 0)) pJ m (£) 


( 2 . 110 ) 
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By employing a Taylor expansion in the neighbourhood of the singular point [68], it can 
be shown 


where 


and 


Also, 


x i~ Vi — (6. ~ Vi) + (£2 - Vi) + 0 (p 2 ) 


= p(^ cos « + ^sin«) + 0(p 2 ) 


= PA,(6) 

(2.111) 

Ai(9) = cos 6 + ^- sin 9^ , i — 1,2,3 

(2.112) 

f 3 ] V2 


A(0) = |x;[A i (0)] 2 | 

(2.113) 

r 2 (p,9) = p 2 A 2 (9) +0 (p 3 ) 

(2.114) 

r 2 “ p 2 A 2 (0) + ° (p) 

(2.115) 

II 

H 

c^. 

-* 1 

c*. 

II 

+ 

0 

XT 

(2.116) 

Jk — n k J — J^(p) + O(p) 

(2.117) 

r,k Jk = O(p) 

(2.118) 

N a = N a (r]) + O(p) 

(2.119) 


Now, it is possible to give explicitly the asymptotic expression of the singular integrand 
function F^ (2.110) as 

F ij (p,9) = -f ij (9)+0(p) (2.120) 

P 

Now, all terms in the traction kernel will be considered to find corresponding fij (9) ex- 
pression required for regularising the integral. 

For expression a {j = ^ (r„- r,j r, n ), 

^ n ‘ j = "(^ +o G))(^ + °w)(^ +o(p) ) ow(iV “ ( ’ ,)+ow) 

= O(p) (2.121) 


hence, it is already regular. 
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For expression ^ , 

= p(^f + o(I)) (®+w)(« , W < OH «(,.)) 

= + (2.122) 

Thus, all relevant terms required for Cauchy Principal value integration have been consid- 
ered to find out corresponding terms required for subtracting out, the singularity. 


Vanishing neighbourhood in the local plane 

The explicit asymptotic expression is now derived for the contour of rr]" in the local pa- 
rameter plane in terms of p and 6. The contour of the neighbourhood r, of radius c on S 
is given by the equation 

e 2 = (x k - y k ) {x k - y k ) ; *, y <•• V (2.123) 

that is, in the {- plane 

e 2 = p z A 1 (0) +• 0 (p 3 ) (2.124) 

Therefore, the contour of a™, image of e™, is given by the equation 

P = a{e t 9) 


= z/3(Q)+ 0(e 2 ) (2.125) 

where ft (6) = 1 /A(6) and it carries all the information required to perform properly the 
limit process in {-plane. 


Final formula for CPV integrals 
The original CPV surface integral is 

ts = I Hij (y, x) N a ({ (*)) dS (x) 

J Sp 

= {iC (». *) (( (x)) dS (x ) | (2.126) 

After the introduction of the parametric representation on each element, the singular in- 
tegral becomes 


I s = lim 

e-t-0 



HiAy,x{t;))N a (o j m ( 0 d{ 


(2.127) 
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and in local polar coordinates, 


I s = lim 

e -+0 



l?m 


(p, 6) dpd(9 


(2.128) 


Adding and subtracting on each element around y, the corresponding asymptotic expres- 
sion (6) /p of the integrand function, the original CPV results in an alternative form 
more suitable for actual computation. 



No limit is no longer present. Because, 


to In £ IE/”’ fv = 0 (2.130) 

Therefore, the original CPV integral on S p can be expressed as a sum of double and one- 
dimensional regular integrals, each one in terms of local coordinates. Each term in the 
final sum refers to one element at the time, thus allowing the computation to be performed 
in an element-by-element fashion. The actual computation is performed by using Gauss- 
Legendre quadrature formulae and dividing the element into triangular subelements. In 
the above formula, p m {9 ) is the value of rho, corresponding to the point at the end of 
zone 1 (section 2.5.1) along the ray pertaining to the Gaussian value of 6 in the element 
under consideration. For the remaining two zones (section 2.5.1) on each ray, ordinary 
Gaussian integration is possible as integral is regular. 


Non-linear transformation 

Above regularization technique consists in subtracting the intrinsic singularity to the sin- 
gular kernel obtaining a regular integral that can be computed via standard quadratures 
and adding afterwards the integral of that intrinsic singularity performed analytically. 
But computation of the first integral on the right hand side of the equation (2.129) is 
poor as the integrand is the difference between two infinities of the same order and coef- 
ficient. So a nonlinear transformation is used in this work having null Jacobian at collo- 
cation point to improve the accuracy of integration. Radial variable p of the polar coor- 
dinates is transformed by the non-linear transformation used by Doblare and Gracia [66]. 
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S' 

P 

Q (Ci (Ci. 6). 6(£i , 6) > ( 3 (£1 1 £ 2 ) ) 
Q' (Ci, C 2 ,0) 

Qo 1 (£ 1 , £ 2 ) > ^2 (£° , £ 2 ) , ^3 (£ 1 , £ 2 ) ) 


integration element 

projection of S on the £i ~ (2 plane 

source point 

point on the element, 

projection of Q on the £1 £2 plane 

the closest point, to P on S. 


This direction-dependent transformation is chosen as it results in distribution of Gauss 
points that is almost a perfect circle about the point of collocation. The transformation 

H e (-1,1) p e (0 ,p m {9)) is 


p(M) = 


B 


r , P {~ M) 


(p + 1)^ + 


L (0) 


B 


16 




-M). 


(//. + l)' 1 (2.131) 


with min ep m {9) > B > 0 in order to be montonic. Here, has been used for B. 

Here, r, p {-1,9) is the same as A{9) computed in equation (2.1 13). 'This transformation is 
seen to improve the accuracy of CPV integrals. 


2.6.1 Quasi-singular integration 

In order to achieve a high level of accuracy in the boundary element, formulation, all inte- 
grals involved must be accurately evaluated. Near-singular integrals need to be calculated 
accurately. Cruse and Aithal [8] proposed an algorithm similar to the singularity subtrac- 
tion technique. This algorithm was used and improved upon by Mi and Aliabadi [69] for 
elastostatic crack problems. In this thesis, the same algorithm has been extended to elas- 
todynamic case. In this method, only the boundary variables are expanded and integration 
is carried out partially in real space and partially in intrinsic space. Accurate integration 
results are obtained. 

The integration coordinate system 

The algorithm is based on the Taylor series expansion of boundary variables. The expansion 
is taken at Q 0 which is the closest point on the element to source point P. The integration 
system is set up on the tangent plane of the element at Q 0 . This tangent plane <1 - £2 and 
the integration coordinate system are shown in figure 2.6. The notation used is as follows: 
The normal vector n of the tangent plane at the closest point Q 0 (£i,£°) is given by 


n = Ax B[\A x B 


(2.132) 
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The terms M,f (j = 1 , 2 ) denote the derivatives of shape functions A/„ with respect to § 
and x? are the values of the geometric coordinates x * at the nodal point . 

The tangent plane is constructed using the normal n and the point. Q () . The integration 
coordinate system takes the integration plane as the coordinate plane' on which Ci , (2 lie. 
The vector e 3 = n is taken as the (3 direction in the integration system and e t = A/|A| 
as the Ci direction. Therefore, the (2 direction is obtained from vector product e 2 = 
e 3 x ex/|e 3 x e x J. The origin of the integration coordinate system, O' is chosen as the 
projection of the source point on Ci — C 2 plane. It is obtained by solving the linear equations: 

(■ P-O') A = 0 (2.136) 

(P-O')-B = (J (2.137) 

(Qo - O') ■ (A x B) = 0 (2.138) 

for unknown coordinates of point O'. The coordinate transformation of a space point X 
in the original system to the integration system X' is given as 

= (Cl, < 2 ,( 3 ) 

= (ex • (X - O') ,e 2 • (X - O') ,e 3 • (X - O')) (2.139) 


The Taylor expansion of the displacement field 

The displacements iq at a time-step N are expressed in terms of the shape functions N a 
as 


Ui(( i,6Hu?JV“(&,6) 

The Taylor expansion of u t about Qo(C?,C 2 °) can be written as 


(2.140) 


T = Ui(Q 0 ) + ^A£i + ^AC 2 + ~ — 
d(i d ( 2 v 2 d( 

= R(O + 0(5 3 ) 
where 


f (AC,)i + (A( ' )(a<2) 1 2 (AG)2 + of 

( 2.11 


<5 = v / ((ACi ) 2 + (AC 2 ) 2 ) 

= V((Ci - Ci(Qo )) 2 + ((2 - C 2 (Qo)) 2 ) (2.142) 

The first order derivatives used above 


can be expressed as 
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and the terms are obtained as 




^2 x -1 

/ dui 



/ aft 


J 

1 du L 


d& 2 



(2.144) 


The terms on the right hand side of above equation can be evaluated explicitly through 


the use of shape functions for geometry Q = CfM“(£i,£ 2 ) and the displacement u x — 
ufN a (£ i,&)- Also, |A = Q (£°, xil) and for each element can be evaluated with 
help of equation (2.139). The second order derivatives can be obtained from 



dCi dCi 



Q ^Cl #C2 
d(i 

dc,2 dC,\ , dCi d& 
9^1 



§£2 §^2. 
d£i 



/ d2u i 
5Ci 2 

5Ci5C2 

<Pu} 
\ Ik? 


\ 


/ d^uj §Uj 3 2 Ci duj d 2 (2 

/ dl? dQ dl? dc 2 


\ 


d 2 m __ duj d 2 Q\ __ duj d 2 C2 
d£id& d£i d£i d& &(>2 ^ 1^2 
d 2 Ui dui d 2 Ci dui d 2 ^2 

S g -acTl W ~*toW 



where §& and can be evaluated from Q = ( fM a (&,&)■ Therefore, && 

are computed by using equation (2.145). 


Semi-analytical integration 


Refer to figure 2.6, the distance between points Q and Q' is denoted by A and the distance 
Q' to Q 0 by <5. There is a relationship A oc 6 2 between these variables, since Q 0 is the 
tangent point of a tangent plane of smooth surface S. The near-singular integral for any 
singular or hypersingular kernel T ki j can be written as 


I 


T kij (P,Q)u k (Q)dS(Q) 


'S' 


' S' 


T kij (P,Q)u k (Q)J(Q')dS'(Q') 


T kij (P,Q)u k (Q)J(Q')-T kij (P,Q')D(u k ) d S'(Q')+ / T kij {P,Q')D(u k )dS'(Q') 

J Js' 


In T la 


(2.146) 


The first integral is evaluated numerically, while the second integral is evaluated semi- 
analytically. The first integral can be written as 

I n = f T kij (P,Q'){D(u k )0(5 2 ) + 0(6 3 )}dS'(Q') (2.147) 

J S' 


T kij (P,Q) = T kij (P,Q'){l + 0(5 2 )} (2.148) 

J(Q') = 1 + 0(S 2 ) (2.149) 

u k = D(u k ) + 0{8 3 ) (2.150) 


since 
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If the integration element is flat, S will coincide with the projection S , and Pkij(P, Q) _ 
T kij (P,Q') and J(Q') = 1. However, this is not the case when the integration element is 

curved. 

Numerical integration for I n is performed in the same manner as explained in sec- 
tion 2.5.1 by subdividing the intrinsic element with closest point as the vertex. And 
Gauss-Legendre quadrature is done along radial and theta directions. 

The second integral I a is evaluated on the (( 1 ,( 2 ) plane after transformation to the 
polar coordinate system. The procedure is similar to that of Cruse and Aithal [8], that is, 
^-integration of the polar coordinate is performed analytically, while ^-integration is done 
numerically. 


2.6.2 Detailed procedure for analytical integration 


Any point Q' on (( 1 ,( 2 ) plane can be represented as (p cos 0, p sin 0). Its corresponding 
coordinates in original or global x a coordinate system can be found by solving the equations 
resulting from (e x - ( x - 0'),e 2 • ( x - 0'),e 3 • ( x - O')) = (p cos 0, p sin 0,0). Original 
oordinates are found as a function of p and 9 as follows: 


D 


x 1 (Q')-x 1 (P) 


Determinant of 



(2.151) 


— [pcos 9 (e 2 2 e 33 - e 32 e 23 ) - psin 9 (ei 2 e 33 - e 32 e 13 )] + aq (O') - Xi(P) 
boi(9) + pb n (9) (2.152) 


In the above, for example, e n ,ei 2 ,e 13 are components of vector e x . Also, x „{()’) and 
x a(P) represent x a coordinates of origin of new coordinate system on the tangent plane 
and source point in the original frame of reference, respectively. Similarly, expressions can 
be found for coordinates x 2 and x 3 . 

Distance from source point to receiver point on projected element on tangent plane is 
r ~ + P 2 )) where e is distance between source point and O'. 

Sine e , r,i= Sl23z*i£l } ^ = 1, 2, 3 

where x^Q’) are coordinates of receiver point on projected element and x ; (P) are coor- 
dinates of source point. Putting values of coordinates of receiving point on the tangent 
plane in expression for r, * yields 


r,i ~ ~ (P<h (9) -l- pbn (9)) 


(2.153) 



2.6 Evaluation of Principal- Value Integral 


49 


where 6oi and b Xl are constants that are function of 9 for the particular element. 


An — Al e 31 + T,2 &Z2 + A 3 e33 

= _ [(601 + blip) e 31 + (602 + b\2p) e$2 + (603 + blip) e 33 ] 

r 

= — f(^ 01 e 31 + 602632 + 603633) + p (611631 + 612632 + 613633)] 

r 

= -(b 0n (9 )+pbi n ( 9 )) ( 2 . 154 ) 

r 

Thus, all terms involved in the kernels can be transformed to original coordinate system 
x a as function of polar coordinates p and 9. 

Similarly, 


dui . , dui A „ .1 d 2 Ui , A „ . 2 5 2 Ui A „ A „ .1 d 2 Ui , A A 


£>(*) = ^(Q°) + ^AC! + ^:AC2 + ^(ACi) 2 + 


5(2 


2 5( 2 


5(i5( 2 AClAC2+ 2M (AC2r 


iV Q uf + ^r(pcos9 - (i(Qo)) + ^(psin0 - ( 2 (^ 0 )) + 2~d^^ pC ° S ^ 


+^~(P C ° S ^ - Ci(Qo))(psin 0 - Ca(Qo)) + ^^-(P sin0 ~ CaCQo))' 
= (A£(6>) + pA“ (0) + p 2 A£(0)K, where a = 1, 8. 


(2.155) 


All derivatives in equation (2.155) are calculated at Q a , the closest point to P on S. Term 
r,i D(uj) can be written as 

r,i D ( Uj ) = ^ (6oi + bup) (Aq + pA° + p 2 A“) uf 

= 1(<& + c%p + 4 it ? + 4if 1 )^ < 2 - 156 ) 

Terms like r n r,j D(uk ) can be written as 

r,i r,j D (uk) = ^(6 0 i + 6i i p)(c^ + 4p + 4p 2 + c? i p 3 )< 

= ^ {d£ (i, j) + d? (i, j) P + d£ (i, j) p 2 + d£ (i, j) p 3 + d“ (i, j) p 4 } < 

Further r 2 = e 2 + p 2 , where e is a constant for a particular element. It is already seen 
in section 2.5.1 that all temporal terms after analytical integration can be reduced to the 
form (a+6r+c/r 2 ). So all terms of displacement and traction kernels, in the expressions to 
be analytically integrated with respect to p, can be expressed completely as a function of 
p and 9. After analytical integration with respect to p, numerical integration with respect 
to 9 is done as follows. 
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(3 


Figure 2.7: Definition for the angular variable integration. 


Theta integration algorithm 

Fig. 2.7 shows an element projected on the tangent plane and the coordinate system for 
semi-analytical integration. Each side of the element is given by three nodes. Any ray to 
the element starting from O' intersects two line segments, unless O' lies on the projected 
element itself. For every 9, p max and p m i n can be found corresponding to the farthest and 
closest point to O' on the element boundary. Now, intersection point of the ray with a side 
of the element is found by solving 


Ci 

= M(0Ci, 

(2.157) 

C 2 

II 

^ry 

to 

(2.158) 

y/x 

= tan# 

(2.159) 


Hence, values of intrinsic coordinate 0 pertaining to points of intersection of ray with both 
line segements of element boundary, are found corresponding to Gauss point value of 9. 
Thereafter, values of £i and £2 can computed from the values of intrinsic coordinates of 
intersection point. And p is calculated by p 2 = + (f. Thus, values of p min and p max are 

computed for any value of 9. 



2.6 Evaluation of Principal- Value Integral 
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Gaussian quadrature formula for integration with respect to 9 is: 

rPmax 1 PPmax 

/ / f(p,9)d9dp = H i -{9,-e l )w i / /(M)dp 

J p min J J Pmin 

= ^q(0 2 -9 1 )w i g(9 i ) (2.160) 

Here, Wi are the weighting factors and are corresponding Gauss points over the interval 
9i — 9 2 . Thus, integration with respect to 9 can be performed numerically. 

This approach for nearly singular integrals is resorted to only if the distance of source 
point to receiver element is less than the length of diagonal of receiver element. 

Angular variable transformation 



Figure 2.8: Nomenclature for angular transformation. 

The number of integration points required in the angular direction to produce accurate 
results using polar coordinates increases rapidly as the source point projection approaches 
the edge of the projected element. This causes a problem not only for nearly singular 
integrals but also for singular integrals when using discontinuous elements, for which case 
the source point can be near the element edges. This phenomenon is termed as ‘edge 
singularity’. 

When O' is very close to an element edge (Fig. 2.8), the value of p m m(0) an d hence the 
integrand with respect to the angular variable in the inte gral ^^ 160) can vary drastically 
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as B varies from d x to 6 Z . Thus a large number of integration points are required to perform 
the angular integration numerically. 

In order to weaken the angular near singularity, Hayami [/0] applied the angular variable 
transformation: 


t(B) = ^log{ 


1 + sin(0 - a) 
1 - sin(# - a) 


} 


(2.161) 


where h is the perpendicular distance from point O' to the element edge and a is the angle 
made by this perpendicular with the Ci axis of the coordinate system (Fig. 2.8). This 
transformation weakens the edge singularity in the angular direction as the Jacobian of 
the transformation is 


df? 1 _ cos(0 - a) 

dt Pmin ifi) h 


(2.162) 


Hence, the left hand side of equation (2.160) can be written as 


rPmax rB 3 Rpmax Rt( 03 ) J/J 

/ / f(p,e)mp = / / / (p, 9(f)) -~<Mp (2.163) 

J Pmin ** @1 J Pmin ^(Ol) 

Now, the integral can be written in the same form as the right hand side of equation (2.160) 
for its semi-analytical integration. This transformation is quite crucial for improving the 
accuracy of integration when source point is near the element edges. 


Radial variable transformation 


Consider a model radial variable integral: 


I 




(2.164) 


where r — vV 2 + e 2 for planar elements. A radial variable transformation may be applied 
to weaken the near singularity due to l/r“. Hayami and Matsumoto [70] showed that the 
T -1 / 5 transformation: 


R(p) — — (p + e) 


■ 1/5 


(2.165) 


gives near optimal results. 

Equation (2.164) is transformed by the radial transformation R(p) as: 


I 


-L 


R{pm&x) pm ^ p 




r n d R 


d R 


(2.166) 



2.7 Quarter-point elements 
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The L x transformation gives: 


p(R) = ( -R)~ 5 -e 

M = 5(iJ) ' 6 

Hence, the integral (2.164) can be expressed as: 

I _ ((-^-5 -e)) m (5(-i7) ~ 6 ) 

A( Pmin ) {((-#)- 5_ e )) 2 + e 2 } 

which, in turn, can be transformed as : 


n/2 


dR 


I 


L 


1 {{-R)-*-e)) m {h{-R)-*) A13 *R 


di?— — da; 


■ {((- B )- 5- e )) 2 + £ 2 } 


/n uj. o , 

_^ 2 + p2\ n / 2 da; 


/(re) da; 


r -i 


where: 


R = 2 [{J? (p max ) -R (Pmin)} 3? 4“ -R (/2 max ) "b R (Pmin)] 


(2.167) 

(2.168) 

(2.169) 


(2.170) 

(2.171) 


Above radial transformation may be used if the radial part of the integral on the projected 
element is numerically evaluated, instead of being analytically evaluated. At times, it may 
not be possible to evaluate integration of the radial part of the integral on the projected 
element analytically as in the case of integration of displacement kernel term for traction- 
singular quarter-point elements because of the presence of a square-root term present in 
the traction variable. 


2.7 Quarter-point elements 

In this crack element, the mid-point nodes of the element sides, normal to the crack front 
for eight-noded quadrilateral, or emanating from the crack-tip for the 6-noded triangle, are 
repositioned at the quarter-point locations on the crack side so that the following relation 
between field variables and distance r exist: 

u(t) — A\ + A%\/r + A$r (2.172) 

t (r) = Ai + A^s/t + A%t (2.173) 

The displacement variation given in (2.172) contains the appropriate square-root of r term 
for modelling the near-front displacement field. However, since traction and displacement 
are independently modelled in the BEM, the quarter-point distortion does not reproduce 
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the correct traction singularity ahead of the crack. The inclusion of the correct singular 
term for the eight-noded element is accomplished by multiplying the shape functions for 
the element tractions by cf>(r ) as given below [71,72], 

^ (r) = rfi (2 - 174) 

The function <f>(r ) is o(^) for r sufficiently small and assumes the value of one at the 
element side opposite to the crack-front or to the crack-tip, for the 8-node or 6-node 
elements, respectively. This ensures that correct singularity is achieved at the crack, and 
traction continuity is preserved along the element boundaries. 

The combination of quarter-point distortion and shape function modification results 
in a quarter-point, traction-singular boundary element characterized by the following dis- 
placement and traction variations 

u(r) = A l +A 2 Vr + A 3 r (2.175) 

t(r) = + (2.176) 

Consider one side of a quadrilateral element, whose mid-side node has been shifted to 



Figure 2.9: Quadratic quarter-point element 

quarter point position (Fig. 2.9). The singularity is included in the representation of the 
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tractions by using modified shape functions: 

= w{-+w{- + w{- (2.177) 
U = t^<j) x + 1%4> 2 + (2.178) 

where 0 1 = |f(£ — 1), <fi 2 — 1 — £ 2 , 0 3 = |£(£ + 1). And and ^ are the modified shape 
functions which include the r -1 / 2 singularity. 

Now, t\ stands for the nodal values of t* divided by the nodal values of </A 


if 

= *f 

(2.179) 

if 

= tf/2 

(2.180) 

if 

= lim t\<F- 

r— >0 z V l 

(2.181) 


The stress intensity factors can be computed directly by the boundary element code as 
follows: 


Kj = t\s/2-Kl 

K n = t[V (2.182) 

Km = 

It may be noted that the use of traction nodal values of the singular element at the crack 
tip for SIF computation is substantially less sensitive to the discretization than any of the 
displacement correlation procedures [60]. 

Almost all details of the formulation have been explained in this chapter. Same formu- 
lation works very well for quasi-static problems as well as by assigning very slow loading 
rates, dynamic problem reduces to static problem. 



Chapter 3 

Quasi-static Fracture Analysis 


The effect of free surfaces on fracture characteristics for elastostatic problems is investigated 
using three-dimensional boundary element method. Stress intensity factors are calculated 
in BEM using traction values at the nodes along the crack front. A close correspondence 
between Kj computed from traction values and Jf integral is seen. Thus, it is possible 
to relate the evaluated state of stress at the crack vertex, where crack front meets the 
free surface, with a value of SIF. Thickness of boundary layer where free surface effects are 
predominant is estimated. Also, it is possible to estimate accurately the critical intersection 
angle of the crack front with the free surface at which square-root singularity is restored 
at the crack vertex. It may become possible to capture free surface retardation of crack 
growth accurately in crack propagation studies with the methodology used in this thesis. 

3.1 Introduction 

A crack front intersecting the free surface in a three-dimensional (3-D) body has been a 
subject of interest for the last three decades. Linear elastic fracture mechanics parameter 
stress intensity factor (SIF) is widely used in the characterization of material properties in 
the presence of a crack. SIF is defined as the strength of the square-root singularity at the 
crack tip. But SIF loses its significance at a corner point where a crack intersects a free 
surface. The power of the crack tip singularity changes in the vicinity of a corner point, 
and it is only possible to define stress intensity factors in an asymptotic sense. Precise 
behaviour is a function of Poisson’s ratio (u) and the inclination of the crack front to the 
free surface [2]. 

Benthem [1] analyzed a problem of a quarter-plane crack in a semi-infinite half-space 
subjected to symmetric loading. He derived that at the corner point, stresses have r -0 - 452 
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singularity for Poisson’s ratio equal to 0.3. Although it does not seem a diamatio change, it 
implies that the stress intensity factor at this point is equal to zero (or does not, exist) since 
it is defined as the strength of the square-root singularity at the crack tip. It also causes the 
energy release rate to be equal to zero, since this can have a bounded, non-zero value only 
in the case of an exact r _1 singularity of the strain energy. Bazant and Estenssoro [2] con- 
firmed Benthem’s values for the singular eigenvalue at the crack vertex through a 3-D finite 
element analysis developed for conical regions. In addition, they considered intersections 
other than normal and found angles at which square-root singularity is restored. They 
used extrapolation procedures to estimate singularity indices corresponding to infinite de- 
grees of freedom from results of eigenvalue analysis of a set of uniform grids. The success 
enjoyed by these local numerical analyses is not found to the same extent in finite element 
treatments of more global or complete problems which contain crack /surface intersections. 
Different predictions were made concerning the variation in energy release rate ((7) with 
some results showing an increase in G in the free surface region, some a decrease and still 
others no change. Burton et al. [4] used a special set of integral equations for an elastic 
half-space weakened by a crack of finite width and infinite depth subjected to transverse 
excitation. They confirmed the existence of a weaker singularity at global crack /surface 
intersections with both integral equation formulation and finite element analysis of a finite 
geometry with adjustable singularity index at the crack vertex. They found by FEM a 
decay of G as free surface is approached (z/a — > 0) with the drop between the maximum 
value and that of the closest Gauss point (at z/a — 0.015) to the free surface being 30% 
of plane-strain value G p when Poisson’s ratio is 0.4. Nakamura and Parks [17] estimated 
thickness of boundary layer region as 3% of plate thickness. They used domain integral 
method to calculate local energy release rate J and calculated a drop in value of J of the 
order of 30% of two-dimensional plane-stress value of J. In FEM, calculation of ./-integral 
at various positions along the crack front is constrained by the fact that, accurate values of 
stresses are available only at Gauss points of the elements. Hence, accurate determination 
of thickness of free surface boundary layer becomes difficult. 

It can be observed in experiments that an initially straight-fronted crack in specimens 
such as a compact-tension specimen usually grows to a slightly curved shape under a fatigue 
load. Three-dimensional LEFM implies that for v ^ 0 the straight profiles associated with 
the varying energy release rates would not propagate as such but rather take up curved 
profiles with constant energy release rates. Burton et al. [4] generated critical crack profiles 
using FEM with each profile having less than 2% numerical variation in G with depth. 
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The crack profiles generated by them were significantly less curved than those observed 
in practice [5]. Similarly, using FEM, Bakker [18] predicted a relative tunnelling depth of 
2.5 — 3% based on constant stress intensity distributions for self-similar crack propagation. 
Experimentally observed tunnelling depths on fatigue cracks are generally larger, closer 
to 5%, which is close to circular front. Using local eigenvalue analysis with FEM, Bazant 
and Estenssoro [2] found the critical intersection angle to be 101°, a value also observed 
experimentally [5] for v = .3. But in global situation, finite element analysis is not able to 
obtain an accurate stress state adjacent to the crack front in the region close to the free 
surface, especially for cracks which are not normal to the surface. Using FEM for analysis, 
Lin and Smith [6] found that slightly non-orthogonal intersection caused a great loss of 
path-independence in the value of J-integral up to 47% near the free surface. 

The present work uses boundary element method to study the effect of free surface on 
state of stress near the crack vertex. SIFs can be computed along the crack front from 
nodal values of traction in boundary element analysis [73]. Amestoy et al. [74] derived for 
three dimensional configurations a path-independent parameter, J”, equal to the potential 
energy release rate per unit of crack extension. The BEM is ideally suited to the evaluation 
of Ji integral, since the required stresses, strains and derivatives of the strains can be 
accurately calculated at the internal points in the body, with respective boundary integral 
equations. So it is possible to calculate J" integral very close to the free surface, helping 
in more accurate estimation of thickness of boundary layer where free surface effects are 
predominant. In the present work, it is seen that SIFs, computed using nodal values of 
tractions at the crack front, closely correspond to J" integral values along the crack front. 
Critical intersection angle of crack front with free surface computed in this work for a 
finite specimen matches closely with the value predicted by Bazant and Estenssoro for a 
half-space problem [2]. 

3.2 The Boundary Element Method 

This section describes the salient features of the BEM applied to the present problem. The 
boundary integral equation relates displacement Ui and tractions ti at the surface S of an 
elastic, homogeneous and isotropic body. It can be written as 

C ij (P)u j (P) + [ Tij ( P , Q) Uj ( Q ) d S= f Uij (P, Q ) tj(Q)dS (3.1) 

Js Js 

where i,j range from 1 to 3 and P and Q are points on the surface S. The free-term 
tensor Cy is calculated by general close-form expression developed by Mantic [56]. The 
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tensors U %J {P,Q) and Ty(P,Q) are the fundamental displacement and traction solutions, 
respectively, for a point force in an infinite domain (i.c., Kelvin s solution). 

In this work, an approach [62] which combines the advantages oi continuous and dis- 
continuous elements in a single formulation has been used. In this iornmlation, expression 
of the shape and interpolation functions are unified in a form suitable lor f ully continuous, 
fully discontinuous or transition elements (i.e., discontinuous along one or more edges). One 
characteristic of the boundary element formulation is that discontinuities of the variables 
across the element edges do not invalidate the convergence of the technique [63]. 

The surface is divided into m eight-noded elements and on each element, the geometry, 
displacement and tractions are represented in terms of the respective nodal values xf, u- 
and ff (s = 1, 2, 3, ... 8) as follows: 

= <!>*(£, v)x* (3.2) 

«<(&»?) = $(£>»/)< (3.3) 

ti(Z,v) = #(&>?)** (3.4) 

where <t > s , $ are shape functions of the intrinsic coordinates £ and ?/. 

The integral equation (3.1) is discretized as 

ra 8 m 8 

+EE =EE <W m 

1 = 1 S =1 i -1 S= 1 

where the surface point P has been taken to coincide with collocation node a of the 
discretization. Coefficients H$ s and Gf/ are given by the integrals 

G » = ’)))«(?,';)■/'(«, (3.6) 

= f_J\up\Q\u))m,’i)-nuiy\^, (a.?) 

where J(£, rj) is the Jacobian. 

Three types of boundary nodes are identified depending on the function these fulfil: 

1. Geometrical nodes: those used for interpolating the boundary geometry. 

2. Interpolation nodes: the points on the boundary which define the interpolation of 
the variables. The system nodal values correspond to these nodes. 

3. Collocation nodes: where the boundary integral equation is enforced. 
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Collocation and interpolation nodes for continuous elements coincide with the geomet- 
rical nodes. Collocation nodes and interpolation nodes for discontinuous elements and 
discontinuous side of transition elements are situated inside the element. For quarter-point 
element, geometrical and interpolation nodes coincide, but collocation corresponding to 
nodes on crack front is done at nodes situated inside the elements on either side of the 
crack front. To reduce the number of equations and enhance accuracy, equations corre- 
sponding to collocation nodes before and after a given interpolation node are added up [64]. 

Writing equation (3.5) for each of the n collocation nodes of the mesh leads, when 
proper boundary conditions are applied and the terms reordered, to a non-symmetric and 
fully populated system of 3 n equations in 3n unknowns. Solving the system of equations 
generated from collocating displacement integral equation at all collocation nodes, bound- 
ary data of displacements and tractions get generated on the full surface of the specimen. 

3.3 BEM Crack Elements 

In this crack element, the mid-point nodes of the element sides, normal to the crack front 
for eight-noded quadrilateral, or emanating from the crack-tip for the 6-noded triangle, are 
repositioned at the quarter-point locations on the crack side so that the following relation 
between field variables and distance r exist: 

u (r) = Ai + A 2 \fr + A z r (3.8) 

t (r) = Ai + A 2 \pr + A$r (3-9) 

The displacement variation given in (3.8) contains the appropriate square-root of r term for 
modelling the near-front displacement field. However, since traction and displacement are 
independently modelled in the BEM, the quarter-point distortion does not reproduce the 
appropriate traction singularity ahead of the crack. The inclusion of the correct singular 
term for the eight-noded element is accomplished by multiplying the shape functions for 
the element tractions by </>(r) as given below [71,72], 

= (3.io) 

The function is o(^) for r sufficiently small and becomes unity at the element side op- 
posite to the crack front or to the crack-tip, for the 8-node or 6-node elements, respectively. 
This ensures that correct singularity is achieved at the crack, and traction continuity is 
preserved along the element boundaries. 
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The combination of quarter-point geometry and shape function modification results 
in a traction-singular boundary element characterized by the following displacement and 
traction variations 

u(r ) = Aj + A%\fr -t- A$r (3.11) 

t (r) = Bx/y/r + B-x + B 3 \fr (3.12) 


Consider one side of a quadrilateral element, whose mid-side nodes have been shifted 
to quarter-point position (Fig. 3.1). The singularity is included in the representation of 
the tractions by using modified shape functions: 




(3.13) 

u = tlft+tfft + iift 


(3.14) 

where ft = |£(£ - 1), ft = 1 - ft, ft = j£(f + I). 



Now, t- stands for the nodal values of tj divided by the nodal va 

lm\s of (f > } : 


? = <? 


(3.15) 

n = tt /2 


(3.16) 

4 = 


(3.17) 


The stress intensity factors may be computed directly by the boundary element code 
as follows: 


Kr = 

K u = (3.18) 

Kin — £3 

It can be noted that the use of traction nodal values of the. singular element at, the crack 
front for SIF computation (3.18) is substantially less sensitive to the discret ization than 
any of the displacement correlation procedures [73], 

Further, for pure mode-I loading of the crack, the local energy release rate G (or J- 

integral) along the crack front can be expressed in terms of local stress intensity factor Ki 
according to 


r _ K](l — ft) 


(3.19) 


Since a plane strain state exists in the limit r -4 0 all along the crack front, the plane 

strain modulus E/( 1 - n ) is used in this expression independent of the far field stress 
state. 
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3.4 Jf - integral 

This section describes a J-integral suitable for 3-D crack configurations. Amestoy et al. [74] 
defined a J" integral for linear or nonlinear elastic materials. Jf integral is equal to the 
potential energy release rate per unit of crack extension, i.e., the negative of the change of 
the potential energy of the configuration when a unit of crack front extends in the direction 
of the principal normal by an amount da at point s only, remaining elsewhere its original 
length. The value of Jf varies with the position s along the crack front and the term Jf(s) 
is used to denote this dependence on position. 

Given a configuration containing a traction-free crack with local axes X\,x 2 and x 3 
attached to the crack front at a distance s along the crack and where x\ and x$ lie in the 
plane of the crack at s, a: 3 is tangent to the crack front at s, and xi is directed into the 
solid (Fig. 3.2), J™(s ) is given by 

Ji(s) = j {W e ni - u it iti) dr - J (ffi3iij > i) t3 cLA 

= Jc(s) + J A {s ) (3.20) 

where 

Jc{s ) = a line integral evaluated over a contour that lies in the principal 

normal plane of the crack front at s and that encloses the crack-tip 
Ja(s) — an integral evaluated over the planar area (surface) enclosed 
by the contour that includes the crack tip 
t{ = components of the traction vector defined according to the 
outward normal n along (7, i.e., ti = a^rij 
Oij = components of the stress tensor 
Ui — displacement in the Xi direction 
W e = elastic strain energy density 

Carpenter et al. [75] and Dodds et al. [76] had demonstrated the path-independence prop- 
erty of J™ integral for 3-point bend specimens. Dodds [77] demonstrated path-independence 
for J” for tension and bending specimens. These studies concluded that quantity Ja is a 
well-behaved function and is less affected by numerical inaccuracies. 

J" is less sensitive to inaccuracy in computations in crack front elements as it is a 
global parameter. Values of stresses, strains and derivatives of strains are computed at 
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interior points using respective boundary integral equations for them. Interior variables 
are calculated more accurately in BEM because these are calculated with the help of 
global equations, rather than by local interpolation as in FEM. Rigby and Aliabadi [78] 
demonstrated the effectiveness of the boundary element method in combination with the 
J-integral technique for three-dimensional problems. 

A typical distribution of Gauss points for Jf calculation is shown in Fig. 3.3. It can be 
seen that the area enclosed by the J-integral contour is divided into a number of segments, 
one of which is highlighted. For area integration, 2x2 Gauss points arc' chosen over each 
segment. J c is evaluated by numerically integrating the contour integral over the outer arc 
( 3 points Gauss quadrature is used for line integrals). The area integrand can be written 
as 


I A = 


doj 3 duj 
dxz 3.x i 


+ C7£3 


0 2 u t 

dxidxs 



I a has a 1/r singularity but the area integral would tend to zero as r 0. 


3.5 Numerical evaluation of boundary integrals 

Analytical integration of the matrix coefficients in boundary integral expiations is not, in 
general, possible and therefore, numerical quadrature must be used. The principal-value 
integrals involved in computation of H matrix (3.7) have been determined with the help 
of algorithm given by Guiggiani and Gigante [7], and non-linear transformation suggested 
by Doblare and Gracia [66]. Algorithm proposed by Cruse and Aithal [8] and further 
extended to hypersingular kernel cases by Mi and Aliabadi [9], has been used for quasi- 
singular integrations. This algorithm is analogous to singularity subtraction technique. 
Nearest point on the element from collocation point is found by Newton-Raphson method. 
The element, over which integration is to be done, is projected on to the plane tangent at 
the closest node. In this method, boundary variables are expanded about, the closest node. 
Integration is partially carried out in real space and partially in intrinsic space. Intrinsic 
space integration is fully numerical. Integration in real space, i.e., on projection of the 
element, is done semi-analytieally in two stages: radial integration is done analytically and 
angular integration is done numerically. Angular variable transformation of Hayami [79] is 

used to weaken the angular near-singularity arising when collocation node is near the edge 
of an element. 
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A new algorithm for G-matrix for traction-singular elements 

In the above quasi-singular integration procedure, radial integration in real space for 
traction-singular, quarter-point elements for G matrix (3.6) is done numerically as ana- 
lytical integration is not possible due to an additional term \fljx in shape functions for 
traction variables. The integrand, having the square root singular term in traction variable, 
is first integrated by parts resulting in two terms: an exact term and another term which is 
numerically integrated using L -1 / 5 transformation given by Hayami and Matsumoto [70]. 
Integration by parts removes square-root singularity of traction variable from the integrand 
as shown below. 

Suppose ABCD is a quarter-point element(Fig. 3.4). Side AB represents the crack 
front. Suppose this side AB is represented by £ = — 1 in intrinsic coordinates. Traction 
variable over the quarter-point element is represented in terms of the nodal values if(s = 
1,2,3, ...8) as 




( 3 . 22 ) 


where N s are shape functions of the intrinsic coordinates £ and 77. And </>(£) = ^ 

Terms l and x are explained in Fig. 3.4. For straight edge elements, x/l can be expressed 
as x/l = a(9)p + b(9 ), where a and b are functions of 9 only. For radial integration in real 
space, a general term'of the integrand having y/ljx term can be expressed as 


I 


/ 


(p 2 + e 2 ) n l 2 yfap + b 


dp 


(3.23) 


A field point on the projection of the element is located by (p, 9) coordinates. Distance 
between field point and collocation node is given by \/ p 2 + e 2 , where e is the distance 
between the collocation point and its projection on to the tangent plane, i.e., center of the 
real space integration coordinate system. 

Integral (3.23) can be integrated by parts and written in the form 


2 y/ap + bp m C 2y/ap + b (p 2 + e 2 ) n/2 mp m 1 - p m n(p 2 + e 2 ) 2 l p , 0 . 

a(p 2 + e 2 )”/ 2 J a (p* + e 2 )" P [ } 


Integral on the right hand side has been numerically integrated using L ~ 1 / 5 transformation 
given by Hayami and Matsumoto [70]. Integration by parts removes square-root singularity 
of traction variable from the integrand, thus improving the accuracy of integration. 
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3.6 Results and Discussion 

Fig. 3.5 shows the geometry of the center crack tension specimen and tin* loading used 
in numerical analysis. The specimen has the parameters: afw O.LM ./»/»■ 1 ,b/u> = 

0.5 where a is the crack half-length, w is specimen half-width, h is half-height, and b is 
specimen half-thickness. The material is assumed to be homogeneous and isotropic with 
shear modulus /x = 90 G Pa and Poisson’s ratio v = 0.3. The center crack specimen is 
loaded by a tensile load at y — dt/i planes. 

Due to the symmetry in the specimen geometry and loading, only one-eighth of the 
specimen is modelled in the boundary element analysis. Fig. 3.6 shows the boundary- 
element idealization of one-eighth of the center-crack specimen. Fight-node quadrilateral 
elements are used. A combination of continuous and discontinuous elements is used to 
reduce the overall degrees of freedom. The eight-node elements bordering the crack on the 
side of the crack surface are quarter-point, elements, while the ones on tin 1 opposite side, i.e., 
ahead of the crack front are quarter-point, traction-singular elements. Triangular elements 
(collapsed eight-noded elements) are used to surround the. intersection point, of the crack 
front with the lateral facets. These triangular elements are modelled as quarter-point, 
traction-singular elements. 

Appropriate boundary conditions are applied to the symmetry planes. 'Phis is done by 
restraining nodal displacements in the direction normal to the plane or planes of symmetry 
to which the node belongs. The nodes on the crack surface, except for those directly on 
the crack front or on the edges with lateral facets, are left free, to move in three orthogonal 
directions. 

Table 3.1 shows details of some of the meshes with varying degrees of refinement. Runs 
were performed for convergence study for a through straight center-crack specimen under 
Mode-I loading. The profiles of K r along the thickness (Fig. 3.7) are seen to converge across 
the thickness except for the free boundary layer region. SIFs are calculated from values of 
traction at the crack front nodes. SIF value at the crack vertex is found to be sensitive to 
the boundary element idealization. The boundary element mesh used for analysis models 
square-root singularity all along the crack front including the crack vertex. But due to the 
actual singularity at the vertex being different from square-root, convergence in SIF at the 
vertex is not obtained with mesh refinement. But an increasing reduction in SIF at the 
vertex is seen, implying that singularity is weaker than square-root. Influence of weaker 
singularity is very well reflected in the lower value of SIF obtained at the crack vertex. 
Value of traction at crack vertex exhibits an integrated influence of reduced singularity 
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over the entire element at the crack vertex. Despite being mesh dependent, value of SIF 
calculated from crack vertex tractions reflects the level of singularity at the vertex. This 
feature may be used to capture accurately the changing singularity trends as intersection 
angle of the crack front changes to non-normal values. 

The 3-D characteristics of the straight crack are studied using the mesh no. 4 (Ta- 
ble 3.1), which has 3317 nodes, 1075 elements and 9951 degrees of freedom. The selected 
mesh has 12 layers of elements on face containing the crack. The thicknesses of different 
layers after normalization with half-crack length are .008, .008, .024, .032, .048, .104, .208, .288, 
.316, .316, .316 and .316 respectively. Various three-dimensional aspects of fracture charac- 
teristics of a through-thickness center crack in a plate are enumerated as follows: 

1. Path- independence of integral: Circular contours around the crack front on 
planes parallel to the free surface are used for J" computations. The computed values 
of Jc, Ja and J" are shown in Fig. 3.8 and Fig. 3.9 as a function of radial distance of 
contour from crack front for several positions on the crack front. Although Ja varies 
with radial distance of outer contour from crack front, r , it clearly approaches zero 
as r -> 0. Also contribution of Ja term reduces sharply on moving away from free 
surface. 

The parameter J™ varies through the thickness of the specimen and is the three- 
dimensional counterpart of the two-dimensional J integral. It should be noted that 
for 3-D configurations, Jc is no longer path-independent as it is for 2-D plane stress 
or plane strain configurations. The term Ja added to Jc gives path- independence 
to J{\ In Fig. 3.9, e.g., for zjb — .004, Jf shows almost zero variation for different 
contours while Jc shows 5% variation. But path-independence of is not preserved 
at positions very close to the free surface. Shih et al. [80] reported a variation of 
25% in the value of point-wise energy release rate near the free surface as calculated 
by domain integral method. In the present study, variations of nearly 50%, 16% 
and 8% are observed in values at z/b = 0.00006, 0.0001 and 0.0002 respectively 
(Fig. 3.8). Here z/b = 0.0 denotes the free surface position. Thus, it is seen that path- 
independence of J 7 / is not preserved as one gets very close to the free surface. At the 
crack vertex, square-root singularity in the stresses is assumed for boundary element 
analysis, while a lower order singularity exists there. Hence, incorrect modelling of 
the stress field might have resulted in loss of path-independence near the free surface 
for Jf integral. But J” values are seen to stabilize for contours at larger radial 
distance from crack front. Remote contour Jf values are taken as correct J" values 
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in this work. 

Fig. 3.10 shows variation of J c along crack front, for cent, ours at. different, radial 
distances from crack front. Towards the mid-plane, J c shows almost no variation 
with radial distance of contour from crack front. Also ./<> matches closely with J 
computed from Kj with equation (3.19) closer to mid-plane. A lot. of variation in 
value of J c , calculated at r/a = .008, is seen in the thickness direction near the free 
surface. It implies stress field is changing sharply in thickness direction near crack 
vertex. 

2. Correspondence of Jf and Kj: It is seen in Fig. 3.11 that ./[' value at the 
free surface approaches the value of J calculated using equation (3.19) from K; 
value computed with traction formula. Thus, SIF calculated from crack-tip tractions 
reflects the state of stress at the crack vertex, though the value of SIF is mesh- 
dependent. 

3. Thickness of free surface boundary layer: Free surface boundary layer is the 
zone adjacent to the free surface where influence of the lower vertex singularity is 
felt. Burton et al. [4] found reduction in energy release rate to be 30% of plane-strain 
value at the closest Gauss point to the free surface (z/a - .015). In results from the 
present study, same reduction in value of J{ 1 is seen at z/a = .001. An FEM analysis 
is constrained by the fact that accurate values of stresses can be calculated only at 
Gauss points of the elements, restricting calculations of G near free surface to the 
plane containing Gauss points of the elements in the layer adjoining the free surface. 
In BEM, values of stresses and derivatives of displacements can be calculated much 
closer to the free surface with the help of the respective boundary integral equations. 
Nakamura and Parks [17] estimated that corner singularity region extended up to 
a spherical radius of about 3% of plate thickness away from the intersection. The 
present study finds that boundary layer thickness, where J\ l varies from low free 
surface values to 2-D plane strain value, is z/a = .008. Thus boundary layer where 
free surface effects are predominant is thinner than its previous estimates. In this 
study, a reduction to the extent of 50% of plane-strain J is observed in J” value close 
to the free surface. Zero value of Jf is not achieved despite a weaker singularity than 
square-root at the vertex because square-root singularity in stresses is assumed all 
along the crack front in the boundary element model used for computation. 

4. Effect of plate thickness on SIF: Fig. 3.12 shows values of normalized Kj along 
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the crack front for v = 0.3 for several thicknesses of the specimens. The two finite 
thickness extremes treated here are: (i) the thick specimen being that used to simu- 
late the basic infinite geometry and having b/a = 2.166, (ii) the thin specimen having 
b/a = 0.091. For both the thick and the thin specimens, value of K[ is lower at free 
surface vertex(z/a = 0.). For the thick specimen, it keeps increasing along the crack 
front till z/a — 0.167 and reaches a peak value there. Thereafter it decreases gradu- 
ally to plane-strain value towards mid-surface. For the thin specimens, intermediate 
peak is never realized, rather Kj keeps increasing in value till mid-thickness. Similar 
behaviour was reported by Burton et al. [4]. 

5. Effect of Poisson’s ratio on SIF: Fig. 3.13 shows influence of Poisson’s ratio v 
on Mode-I stress intensity distribution along the crack front. At the crack vertex, 
Kj reduces with increase in u, implying singularity of stress field at the vertex grows 
weaker with increase in v. Benthem [1] also found that Mode-I singularity exponent 
at the crack vertex gets weaker with increase in Poisson’s ratio. 

Opposite trend is observed with respect to Mode-II loading for normal intersec- 
tion(Fig. 3.14). Ku at the crack vertex keeps increasing with increase in Poisson’s 
ratio, implying an increasing trend in singularity exponent. Oscillation is noticed in 
Ku curves near the free surface. Explanation for this is that for Mode-II loading, 
singularity is stronger than square-root at crack vertex resulting in computed Ku 
value at the crack vertex being higher than that at the mid-surface. Consequently, 
free surface layer in pure Mode-II loading is taking more load or providing more 
resistance to the applied load. This results in inner layers taking less load in order 
to maintain equilibirium. Opposite behaviour is noticed for Mode-I loading where 
inner layers have higher SIF to compensate for part of the crack resistance lost in 
the breakthrough. Perhaps it is more difficult to handle a singularity stronger than 
square-root compared to a singularity weaker than square-root. A still thinner layer 
at the free surface would have smoothened the curve for Ku distribution along the 
crack front. Oscillations are also noticed in Ki curves in Fig. 3.7 for coarser meshes 
and Ki curves are smoothened on further refinement of the mesh. It seems that for 
Mode-II loading, greater mesh refinement is required compared to Mode-I loading. 

6. Transition from plane strain to plane stress in vicinity of crack front: 
Plane strain constraint, defined as o zz jv{a x x + cr yy ), is equal to 0 for plane stress 
and to 1 for plane strain. For thick ( b/a = 2.166) and thin ( b/a = .091) plates, 
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contours in Fig. 3.15 and Fig. 3.16 respectively show how constraint values are lower 
in free boundary region compared to inner layers. Constraint, keeps increasing with 
decreasing radial distance. For the thick plate, ze.ro constraint, legion exists beyond 
radial distance from crack front, r/a = .002 near the free surface', and beyond r/a = 
1.7 near mid-surface. If constraint > 0.9 is used as the criterion for the plane strain 
condition, the boundary of plane strain region lies near a radial distance of .03a from 
the crack front at the mid-plane and at r/a = .0002 for z/b = .003. Near free surface, 
at zjb — .00003, constraint value of .15 is computed at r/a = .0002. 

7. Effects of crack front tunnelling: Through-thickness fatigue cracks in plate 
geometries almost always show crack front tunnelling, i.e., a curved crack front with 
the the point of maximum length of the crack at the center. It is also observed that 
fatigue crack fronts penetrate the free surface at a non-zero angle with the normal 
to the free surface. Bazant and Estenssoro [2] argued that, this angle is fixed by 
the angle at which the square-root singularity of stresses and strains at the free 
surface penetration is restored. They showed that the front edge of a planar Mode-I 
crack that propagates in an elastic isotropic plate in a plane normal to the plate 
surfaces must terminate at the surfaces with a certain angle ft depending only on 
Poisson’s ratio v. Tests were carried out by Rolfe and Barsom [5] on structural steel 
in static loading according to E-399 ASTM Standard. Their test results showed that 
theoretical lines of inclination ft = 101° with the plate surface agreed very well with 
the observed terminal directions of the first arrest marks, which corresponded to an 
essentially elastic fracture stage. 

In the present study, calculations on tunnelling crack fronts are performed using cir- 
cular crack fronts with different penetration angles, viz., 95°, 98”, 101”, 105" and 1 10°. 
The model selected for studying crack-tunnelling effect, has 18 layers of elements of 
unequal thicknesses on face containing the crack. The model has 4245 nodes and 1345 
elements and 12735 degrees of freedom. The thicknesses of different, layers after nor- 
malization with half-crack length are .008, .008, .016, .024, .024, .032, .048, .048, .048, .048, 
.064, .064, .064, .064, .08, .08, .08 and .08 respectively. The specimen has geometrical pa- 
rameters. a/w — 0.24, h/w = l,b/w = 0.22. Also in the present, study, mesh used 
for analysis of curved crack fronts is not orthogonal . For an orthogonal mesh, the 
elements surrounding the crack front have sides normal and parallel to the crack 
front. In a non-orthogonal mesh, the orientation of the element sides is arbitrary. In 
this study, the element sides that would be normal to the crack front for an orthog- 
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onal mesh are taken to be parallel to the free surface. This makes the direction of 
self-similar crack growth as parallel to the free surface in this study. The crack can 
grow in a self-similar manner in two directions: normal to the crack front or parallel 
to the free surface. A crack advance associated with energy release rate G calculated 
in the plane parallel to the free surface produces a crack that retains its shape [4]. 
Accordingly, self-similar crack growth direction has been taken as parallel to the free 
surface in this study. 

Fig. 3.17 shows a comparison of the calculated through-thickness ^/-distributions for 
different penetration angles. Values of Kj are calculated from crack front tractions. 
These values exhibit certain oscillations along the crack front, as the modelled crack 
front is not smooth but is made of small straight line segments. The position of the 
corner nodes of crack front elements is determined from the assumed circular shape 
of the crack front, and the position of the mid-side nodes are then obtained as the 
mid-points of both neighbouring corner nodes. For correct interpolation of traction 
singularity in crack front elements in BEM, a restriction is imposed on the positions 
of nodes in the elements [81] requiring that crack front be straight. SIF values at 
the mid-side nodes is about 10% higher than those at the corner nodes. This SIF 
variation is caused by the inaccurate positioning of the mid-side nodes. The large 
difference of SIF values between the corner and the mid-side nodes demonstrates that 
the SIF variation along the crack front is very sensitive to the crack shape, even a 
small difference in crack front shape will cause a large discrepancy of SIF. Similar 
behaviour was reported by Lin and Smith [6] for a penny crack when polygonal line 
approximation is used for the crack front. They reported an 11% higher value at mid- 
side nodes compared to corner nodes, citing inaccurate positioning of the mid-side 
nodes as a reason for variation in SIF. But their SIF values were estimated by the 
J-integral method. They found it difficult to obtain an accurate stress state adjacent 
to the crack front in the region close to the free surface, especially for cracks which 
are not normal to the surface. They reported a loss of path-independence up to 47% 
near the free surface even for a slightly non-orthogonal intersection. 

Fig. 3.18 shows the variation of Kj for different crack front intersection angles with 
free surface, after ignoring the mid-side node SIF values resulting in smooth profiles. 
It is seen that free surface SIF is greater than mid-surface SIF for intersection angles 
of 101° and higher. For intersection angle of 105° and beyond, SIF at the crack vertex 
is higher than SIF values of inner layers. 
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3.7 Conclusions 

An elastostatic three-dimensional boundary element, analysis of the center-crack tension 
specimen is performed. Values of J” integral calculated along the ci.uk bout correspond 
well with the Ki values calculated from nodal traction values. At. the earner intersection 
point of the crack front with the free surface, values of K h directly available from BEM 
solution, correspond closely with Jf value calculated on plane as close as z/a = .000125 
to the free surface. Increase in Poisson’s ratio reduces K r value at the vertex for Mode- 
I problem and increases K n for Mode-II problem. This trend is similar to Benthem’s 
results [82] where stress singularity became weaker for Mode-I and stronger for Mode-II 
with increase in Poisson’s ratio. 

For curved crack fronts with non-normal intersection with free surface, }(, values at 
crack vertex are higher than mid-surface Kj values for ('.rack intersection angle of 101° and 
higher than those of inner layers for intersection angles of 105" and above. Bay, ant, and 
Estenssoro [2] had also estimated critical angle to be 10,1° for v ■ 0.3 using a finite element 
technique. 

The value of SIF at the crack vertex, computed from crack- vertex tract ion, is seen to be 
mesh dependent. But SIF value reflects the influence of the exponent, of singularity existing 
at the crack vertex. A lower value of SIF is computed for a weaker singularity at. the vertex 
and a higher value of SIF is computed for a stronger singularity at, the vertex. For the 
same level of mesh refinement, value of SIF is seen to be sensitive, to change in singularity 
of stresses at the vertex brought about by changes in either Poisson’s ratio of the material 
or intersection angle of the crack front with the free surface. It is possible t.o estimate the 
intersection angle at which the SIF at the crack vertex is equal to or higher than SIF of 
inner layers. So the intersection angle at which free surface vortex point becomes critical 
along with other parts of the crack front is estimated. 

Burton et al. [4] surmised that the underlying theoretical assumptions and developments 
in LEFM were to a degree incapable of explaining all of the various phenomena actually 
present in retardation of crack front near the free surface. Bakker [18] mentioned the need 
of invoking plasticity and crack closure effects to explain free surface retardation. On the 
other hand, for quasi-brittle behaviour with only small-scale yielding, Smith (19] found 
that cracks identical in shape to those in fatigued metal models could be grown under 
monotonic loads (R = ratio of minimum Kj to maximum Kj = 1.0) in photodastic bodies 
suggesting that closure effects might not strongly affect crack shapes. 

The present study shows that the critical angle of intersection predicted by Bazant and 
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Estenssoro [2] for a half-space problem is also achievable by analysis of finite specimens 
using BEM. So this study removes the doubt that LEFM cannot explain 101° angle as 
critical angle of intersection for a Mode-I finite specimen. The presented approach may help 
in predicting crack growth profiles accurately for problems where crack front penetrates 
free surface. 
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Table 3.1: Mesh characteristics for Mode-1 run s with straight, crack front,. 


Mesh no. 

crack front 
element 
size 
(A d/a) 

free surface 
element 
size 
(A z/a) 

size of regular element, s 
(size/clear distance to 
crack front) 

(on face having crack) 

No. of 
elements 

No. of 
nodes 

1 

0.1 

0.1 

0.5 

•108 

1250 

2 

0.1 

0.008 

0.5 

020 

1900 

3 

0.05 

0.048 

0.25 

1182 

3042 

4 

0.05 

0.008 

0.25 

1075 

3317 

5 

0.05 

0.0008 

0.25 

1239 

3821 
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Area segment 
• Gauss-points 


Figure 3.3: Distribution of internal points about the crack plane for Jf calculation. 
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Figure 3.4: Quarter-point element projected in the plane tangent, to the nearest node. 
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Figure 3.6: Boundary Element Mesh. 




Figure 3.7: SIF distribution along crack front for different meshes (Table 3.1) 
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Chapter 4 

Dynamic Analysis: Stationary Crack 


This chapter deals with the implementation of a three-dimensional time-domain bound- 
ary integral formulation for a center-crack, finite solid under symmetrically applied step 
loading. The BEM displacement time domain formulations have, hitherto, been limited to 
analyzing two-dimensional crack problems, though hypersingular formulations have been 
used to analyze finite cracks in infinite domains. In this chapter, variation of dynamic stress 
intensity factor (DSIF) along the crack front for a stationary, through-thickness straight 
crack is studied for a finite solid under step loading. The state of stress is evaluated at the 
crack vertex, where crack front meets the free surface. The effect of free surface on DSIF is 
investigated. The effect of waves travelling in thickness direction is explained. It is possible 
to estimate accurately the critical intersection angle of the crack front with the free surface 
at which square-root singularity is restored at the crack vertex under step loading. A new 
partitioning scheme is proposed for spatial integration of elastodynamic kernels. 

4.1 Introduction 

While 2-D time domain boundary element formulations have received a lot of attention 
in recent years, 3-D time domain boundary element formulations have seen only limited 
applications. The reason is that time-domain boundary element algorithms become expen- 
sive for calculating response over large periods of time as information from earlier times 
must be used. The computational time and memory requirement become very large for 
3-D problems. 

Two different time-domain integral equation formulations are popular for the applica- 
tion of BEM for fracture problems. One is the conventional formulation where the integral 
representation of the boundary displacements is discretized. The second one is the hyper- 
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singular formulation where integral representations for both displacements and tractions 
are discretized. Hypersingular formulations have been used to analyze finite plane cracks in 
3-D infinite domains [83,84]. Laplace transform BHM has also been used to analyze cracks 
in finite solids under dynamic loading [85]. Wen, Aliabadi and Hooke j8ti,87| presented 
a Laplace transform displacement discontinuity and fictitious stress formulations for 2-D 
and 3-D problems. They determined the stress intensity factors for several 2-D and 3-D 
problems using an equivalent stress approach. Dominguez and Gallego [tip] used the di- 
rect time domain boundary element formulation in combination with the traction-singular 
quarter-point elements that had been used by Martinez and Dominguez [73] and Chirino 
and Dominguez [88] for static and time-harmonic computations, respect ively. Time domain 
formulation is advantageous where an accurate solution is required from the very beginning 
and is the best possible approach for problems with time-dependent geometry [50], For 
this study, direct time-domain approach is chosen so that the approach used to study free 
surface effect for a stationary crack can be easily extended to study free surface effect for 
a running crack in a finite solid. There are certain disadvantages associated with other 
approaches. In frequency domain analysis, internal damping reduces the maximum value 
reached by the SIF as a small amount of internal damping is required in frequency domain 
analysis to avoid spurious oscillations. In the Laplaoo transform approach, the solution 
y to be obtained for a certain minimum number of Laplace transform parameter values 
» order to obtain the inverted solution in the retd time tin. I„ w „f moving 

will a ^ °^ emS 35 111 crack Propagation, calculations lor a set of transform parameters 

til? l! perf0rmed f0r each b0Undary confl e ural »"> thus, increasing the computa- 

remol °“ Sh thC dUai reCIpr0city a PP roac b requires loss computer time, it docs not 
duce very well the quick changes that SIF may have along lime [S>|. 

111?” 1 W ° rk ’ ^ Sta " dard 3 ‘ D integral formulation is 

ofDSIF is st ?a center ' crack ’ fimte solid ,mdor symmetrically applied loading. Variation 
Partitioning h ® ** *• th ™>gl>-thickuc«s crack. In this work, a new 
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show crack front tunnelling, i.e., a curved crack front with the deepest point of the crack 
at the center. Benthem [1] analyzed a problem of a quarter-plane crack in a semi-infinite 
half-space subjected to symmetric quasi-static loading. For a crack front normal to the free 
surface, he derived that at the corner point, stresses have r - - 452 singularity for Poisson’s 
ratio (y) of 0.3. By a local finite element analysis, Bazant and Estenssoro [2] estimated the 
critical intersection angle to be 101° (for u — 0.3) for a stationary crack under elastostatic 
loading. By a finite element formulation similar to that used by Bazant and Estenssoro [2], 
Gudmundson and Ostulund [3] found that a weaker singularity also exists at the corner 
point for a dynamically growing crack and estimated the critical angle to be 101°, almost 
independent of the crack-tip velocity. The success enjoyed by these local numerical analy- 
ses is not found to the same extent in finite element treatments of more global or complete 
problems which contain crack/surface intersections. Different predictions were made con- 
cerning the variation in energy release rate ( G ) with some results showing an increase in 
G in the free surface region, some a decrease and still others no change. Burton et al. [4] 
generated critical crack profiles using FEM with each profile having less than 2% numerical 
variation in G with depth. The crack profiles generated by them were significantly less 
curved than those observed in practice [5]. Thus in global situation, finite element analysis 
is not able to obtain an accurate stress state adjacent to the crack tip in the region close 
to the free surface, especially for cracks which are not normal to the surface. Using FEM 
for analysis, Lin and Smith [6] found that slightly non-orthogonal intersection caused a 
great loss of path-independence in the value of ./-integral up to 47% near the free surface. 
Agrawal and Kishore [90] demonstrated the effectiveness of BEM in evaluating the state of 
stress at crack vertex and estimating critical intersection angle accurately for quasi-static 
loading. They achieved a close correspondence all along the crack front between Jj 1 and 
Ki computed from the tractions at the nodes on the crack front. In a review article, Nish- 
ioka [91] has reported from an FEM analysis that value of the dynamic /-integral reduces 
near the free surfaces of the plate for a through-thickness straight crack. But author is 
not aware of any work in which critical intersection angle has been estimated for impact 
loading for a finite solid. 

The present work uses boundary element method to study the effect of free surface on 
stress intensity factor distribution along the crack front for a stationary through-thickness 
crack in a finite solid under step loading. SIFs can be computed along the crack front from 
nodal values of traction in boundary element analysis [60]. Straight as well as curved crack 
profiles are studied. Critical intersection angle of crack front with free surface computed 
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this work for Mode-I step loading for a finite solid matches closely with the value of 
critical angle predicted by Bazant and Estenssoro |2| for clast ost at ie loading, for a half-space 
problem. 


4.2 Transient Boundary Integral Formulat ion 


Consider a homogeneous, isotropic and linearly elastic body of volume V bounded by a 
regular surface S, with its mass density p and dilatational and shear wave speeds as c x and 
c 2 respectively. Using Stoke’s state of quiescent past corresponding to a point load with a 
Dirac delta function variation, Love’s integral representation can be derived using Betti’s 
reciprocal theorem. For zero initial conditions and zero body forces, t he displacement of a 
point y at time T can be represented by the following boundary integral equation: 


Cij{y)ui{y,T ) 





fI v (x,y,T) * tt,(x. /'!] dS{x) (4.1) 


where 


G%j * t{ 
a%j ^ 


Gij(x,T;y,r) U{x,r) d? 


I Hij(x,T]y,r) •«*(*, 
Jo 


r) dr 


(4.2) 

(4.3) 


are Riemann convolution integrals. In the above, Ui{x, r) and /,( z, r ) are the displacement 
and traction components respectively at point x and time r. 'Hie free-term tensor ty is 
calculated by general close-form expression developed by Mantle (56], The fundtunental 
solutions Gij and H tj are displacements and tractions in direction i at a point x at time 
T due to a unit impulse vector in direction j applied at a point y at a preceding t ime r . 
The fundamental solutions are given in explicit form as 


Gij(x,T -,y,r) 
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1 / 0 ! 
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l/r, 
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is the Dirac delta function. 
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- fij( l ~ 2c2 /cf){8(v -r/ci) + (r/ci)S(v - r/ci)} - gij{5(v - r/c 2 ) 
+(r/c 2 )S(v - r/c 2 )}l 


(4.5) 


where v = T - r; hi = - Vi - d tJ = h i h 3 h m n m /r 5 - = h^/r 3 ; gij = + 

fiijh m n m )/r 3 ; + gij . Also n, are the components of outward pointing normal 

vector and n m /i m = n x /ii + n 2 h 2 + n 3 /i 3 . 


4.3 Numerical implementation 

The time-domain BEM in three-dimensions consists of two basic steps as explained in 
[49,58,59]: 

1 . A discretization of the real time axis into a sequence of equally spaced time intervals 
with the assumption of linear variation of displacement and constant variation of 
tractions over each time interval. When piecewise linear variation is used for both 
displacements and tractions, the time-stepping scheme becomes unstable [60,61]. 

2. A discretization of boundary into 8-noded elements. 

The time span of interest is divided into N equal time steps of size AT so that 

T n = n AT, n= 1,2,..., AT. (4.6) 

On the basis of these discretizations a time stepping solution of equation (4.1) can be 
established for the boundary displacements and tractions over each element and for each 
time step. The displacement boundary integral equation can be written as 

r fTig-i r 

CiMy , T n ) - / / [Gijk - HijUi] dSdr = / / [G^ - H ijUi ] dSdr (4.7) 

J T/v-i J S Jo Js 

where the integral on the right hand side is the contribution due to past dynamic history. 
This equation is also an exact statement since no approximation has yet been introduced. 
However, in order to solve the above equation, one has to approximate the time variation 
of the field quantities, in addition to the usual approximation of spatial variation. The 
time integration in the above equation is done analytically, and the spatial integration is 
done numerically. A system of algebraic equations is obtained at time step N through 
nodal collocation, i.e., by allowing y to coincide sequentially with all the nodal points of 
the boundary. 
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4.3.1 Temporal and spatial interpolation of displacement and 
traction 

The boundary variables are interpolated in terms of their values at nodes by 


N 


Ui(x,T) = £[«?*? V) I */?• 

N 

ti(x,r) = £^C(.) I rt'Cml 


(«) 

(4.S) 


71=1 


and 


u 


,,n— 1 


(x) = j]jv'K(x)]«r, 

C=1 

(X) = £]jvK(*)|«; n 


lie 


i - i.2,d 


/ I, Ad 


(4.10) 

(4.11) 




where u" c and u^ 1 ’ 0 represent nodal values of displacement at time nodes n and (n-1) 
respectively. Similar interpolation is done for traction variables also. Here. .V is the total 
number of time steps. Mj and Mp are the temporal interpolation functions, related to 
local time nodes I (initial) and F (final), and are of the following form for piece-wise linear 
variation: 


where 


or 


M < = lj w- « T > 

M i = Mr) 


<f>n{ T ) 


1 if (n - l)AT • / « m.Y/' 
0 otherwise 


(4.12) 

(4.13) 


^n(r) = H (r - (n - 1)AT) - //(? nAT) (4.14) 

H being the Heaviside function. 

The displacements are interpolated by piecewise linear functions. Tractions are repre- 
sented by piecewise constant functions which allow for sudden jump of their values from 
step to step. For tractions, temporal interpolation functions are of the form 


Mi = 0.5 <t> n (r) 
= 0.5 <j> n {r) 


(4.15) 

(4.16) 
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For spatial interpolation of boundary variables, an approach [62] which combines the 
advantages of continuous and discontinuous elements in a single formulation is used in 
this work. In this formulation, expressions of the shape and interpolation functions are 
unified in a form suitable for fully continuous, fully discontinuous or transition elements 
(i.e., discontinuous only along one or more edges). Three types of boundary nodes are 
identified depending on the function they fulfil: 

1. Geometrical nodes: those used for interpolating the boundary geometry. 

2. Interpolation nodes: the points on the boundary which define the interpolation of 
the variables. The system nodal values correspond to these nodes. 

3. Collocation nodes: where the boundary integral equation is enforced. 

One property of the boundary element formulations is that the discontinuities of the vari- 
able across the element edges do not invalidate the convergence of the technique [63]. 
Collocation and interpolation nodes for continuous elements coincide with the geometrical 
nodes. Collocation nodes and interpolation nodes for discontinuous elements and discontin- 
uous side of transition elements are situated inside the element. For quarter-point element, 
geometrical and interpolation nodes coincide, but collocation corresponding to nodes on 
crack front is done at nodes situated inside the elements on either side of the crack front. 
To reduce the number of equations and enhance accuracy, equations corresponding to 
collocation nodes before and after a given interpolation node are added up [64]. 


4.3.2 Analytical integration of temporal terms 


The time integration in equation (4.7) after using equations (4.8) and (4.9) is done analyt- 
ically and the spatial integration is performed numerically. Explicit expressions for time 
integrals may be found in [49,58]. Temporal integration boils down to integration of terms 
like f AT 5 (T — t — r/c ) Mj (r) dr, f AT A 5(T-r- A r) M} (r) dAdr etc. 

For example, 



— t — r/c) (1— r/AT)0i(r) (r) dr = 


( 1 - 


AT' T cAT 1 


0 


if c(T - A T)<r<c: 
otherwise 

(4.17) 


Similarly, all temporal terms can be expressed in the form a + br + -^ after analytical 
integration in time. Also a, b, d are constant in the separate zones defined in section 4.3.3. 
It is noticed that for the distance reached by both shear and compression waves travelling 
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from a source node up to a time AT, coefficient of r \ is mo. Thus, the temporal terms 
are not adding to the singularity in r. Here, c is compression (e,) or shear (c,) wave speed. 


4.3.3 Details of spatial integration 

The kernels in elastodynamic case have discontinuous derivatives at points such that 
r/d( 2 ) = /cAT (k integer). Therefore, the standard numerical methods, like Gaussian 
integration, give poor results if applied over the whole element., say, Instead one has to 
partition each element Tj into sub-elements : A, = U/v* with /vf - {{k 1 )A/’ < r/c < 

A; AT}. Then, an integral over E{ is the sum of sub-integrals man- Tf, each sub-integral 
being evaluated using standard numerical methods. This technical difficulty is common to 
all BEM development using retarded potentials and time-marching schemes. It is reason- 
ably tractable for 2-D problems. The partitioning is dearly a purely geometrical problem 
(find the intersection points of a given curve C with circles centered at a given pointy 
and of equally spaced radii). On the contrary, its three-dimensional equivalent (find the 
intersection curves of a given surface S with spheres centered at a given point y and of 
equally spaced radii) is extremely complicated. However, Karabalix [;d)| had performed the 
spatial integrations by partitioning the receiver elements into very small rectangular sub- 
elements. With his scheme, the possibility of only certain portions of the source element 
being ‘active’ also could be taken into consideration. 

In this work, a new scheme is proposed for partitioning of surface elements for 3-D 
problems as follows. 


Partitioning scheme 

Consider a boundary element S p . It is mapped onto a square region H (Pig. 4.1) in the 
parameter plane having local intrinsic coordinates £ » (d ,(■>)■ Suppose source point is 

y and a typical receiver point is an Hence, a typical integral over the element S p can be 
written as 


Is i (*(£)» r )N c (£, as) M(r)dr <\S(x) 


i'T (£,*) d.V(x) 


= f \f {*{£), T-,y,T)M(r)dT 

J Sp uo 

= LI [Vij ( ' X ^' T 5 y ’ r) M ( r ) dr l N " (0 J (0 dCdC- (4-18) 
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where Vij is displacement or traction kernel, M(r) is temporal interpolation function, and 
iV c (£, x ) are spatial interpolation functions of field variables. 

Coordinates of a receiver point on a boundary element are given by a parametric rep- 
resentation in terms of shape function M c (£ 1; £ 2 ) and nodal coordinates x\. 


x i = Y J M c {^2)x c i , i = 1,2,3 (4.19) 

C 

Let us indicate by r] = (771 , 772) , the image of the point z on the element. Point z is closest 
to the source point. The first, second and third order derivatives are obtained by 


dxi dM c 


dfik d£ k 


■Xt 


d 2 Xi d 2 M c 


d 3 Xi 


d z M' 


-X; 


-X- 


(4.20) 


(4.21) 


(4.22) 


By employing a Taylor expansion about the closest point, it is easy to establish formula of 
the form 



dxi 


v dXn 
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3 £i< 9£ 2 
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(£1 - - m) + d g 


(£2 - ihY 




+ J(£i d Xl 




2 V ^ ,i/v ^ 


Z=r) 


(4.23) 


Since geometry is interpolated with help of quadratic functions, above formula is exact. 
The terms and do not appear in the above expression as these are equal to zero 

Cm 

due to quadratic interpolation of the geometry. 

Following a common practice in the BEM, polar coordinates (p, 9 ) centered at rj are 
defined in the parameter plane as 


£1 = 7?1 + p COS 6 

£2 = V2+P sin 9 


so that d£i d £ 2 = p dp d 9 . Hence, 


= e* + p 


dxi 

Wi 


dxi 

W2 


cos 9 + sin 9 


(4.24) 

(4.25) 


Xi - Vi 
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r 2 = E(xi-pi ) 2 

= 2e? + 2 pEiMW + P 2 (ZA 2 (0) + 2B(c,B i (9))) t <- 

+p 4 (22(4W(^)) + ££?(*)) + 2p*Z(B i {0)C i m t />* ?(*) (4.2?) 

For integration, the integration square (-1 < £, < 1) is subdivided into t riangles with 
common vertex at rj = (%,%)• The integral (4.18) of a. typical kernel over each triangle is 
calculated as 

/'(h I'fhnux 

1 = 1 I" (4.2!) 

where 

Fij (p, 9) = Wn (y, x (£ {p, 0))) N r ($ (p. 0)) /»./({ (p, P)) (4.29) 

where Wy = f Q Vij (x,T;y,r) M(r)dr. The limits 0\ and ()■> art* angles subtended by the 
arms of the triangle on the vertex V = (, lum ). Here p nw , is the value of p on a particular 
ray making angle 9 , where it meets the opposite side of the triangle. 

First of all, Gauss integration is performed over 6 direction. Corresponding to the value 
of 6 computed as per Gauss points, p max is calculated for the particular ray. The values 
of p on a particular ray corresponding to the distance of rjnAT or are computed 

by a combination of bisection and Newton- Raphson methods applied on equation (4.27). 
For each ray, the value of p are found corresponding to positions of heads and tails of the 
wavefronts of compression and shear waves, starting from the source nolle, for the element 
and time interval under consideration. The ray is divided into three different zone, for 
separate Gaussian integration. The three different zones are: 

1. both compression and shear wave cover part of the ray. 

2' onl y compression wave covers part of the ray. 


• only shear wave travels that portion ot the ray in that particular tin,,, interval. 

In a time interval, only those areas contribute to the integral which were not affected 

such Z W Z m PrCTi0US time intemk Tlms . jump in the kernels at points 

such that r = c m nhT is handled accurately. 
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Evaluation of Principal- Value Integral 

Principal value singularity occurs for traction kernel Hij when source element becomes also 
a receiver. All temporal terms can be expressed in the form a + br + £ after analytical 
integration in time, as explained in section 4.3.2. Also a, b, d are constant in the respec- 
tive zones defined in section 4.3.3, allowing use of techniques developed for integrating 
elastostatic kernels by splitting the radial integration over three separate zones. 

The Cauchy principal value integral is treated in the same way, as done by Guiggiani 
and Gigante [7] for the elastostatic case. The non-linear transformation suggested by 
Doblare and Gracia [66] is also incorporated. 

Quasi-singular integration 

Nearly-singular integrals need to be evaluated accurately. Algorithm proposed by Cruse 
and Aithal [8] and further extended to hypersingular kernel cases by Mi and Aliabadi [9] is 
applied to elastodynamic case. This algorithm is analogous to singularity subtraction tech- 
nique. Nearest point on the element from collocation point is found by Newton-Raphson 
method. The element, over which integration is to be done, is projected on to the plane 
tangent at the closest node. In this method, boundary variables are expanded about the 
closest node. Integration is partially carried out in real space and partially in intrinsic 
space. Intrinsic space integration is fully numerical, while integration in real space, i.e., on 
projection of the element, is done semi-analytically in two stages: radial integration is done 
analytically and angular integration is done numerically. Angular variable transformation 
of Hayami [79] is used to weaken the angular near-singularity arising when collocation node 
is near the edge of an element. For traction-singular elements, integration of G matrix is 
performed using the procedure explained in [90]. 

4.3.4 Crack Elements 

Quarter-point elements are used along the crack front and traction-singularity is incorpo- 
rated in the quarter-point elements ahead of the crack front [71,72]. The stress intensity 
factors can be computed directly by the boundary element' code as follows: 

Ki = f 2 V^d 
Kji = i\ V 27 7l 
Kni = il\/2nl 


(4.30) 
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where t{ are the tractions (with square-root singularity embedded in them) ah nodes on 
the crack front and l is the length of crack front elements in the direct ion of singularity, 

It may be noted that the use of traction nodal values of the singular element, ah the crack 
tip for SIF computation is substantially less sensitive to the discretization than any of the 
displacement correlation procedures [60]. 

4.4 Results and Discussion 

Elastic analysis of a stationary crack subjected to symmetrically applied step load has 
been performed for this work. Fig. 4.2 shows the geometry of the center ('rack tension 
specimen and the loading used in numerical analysis. The specimen has the parameters: 
a = 0.24 cm, to = 1 cm, h = 1 cm and b = 0.5 cm where a is the crack half-length, to 
is specimen half-width, h is half-height and b is specimen half-thickness. The material is 
assumed to be linear elastic with properties ms: shear modulus // - 90 ( II ’a, Poisson’s ratio 
v = 0.3 and density p = 7800 kg/m 3 . Hence, c.\ = 6354 m/s and <■■> 3390 m/s. The 

center crack specimen is loaded by symmetrically applied tensile step pulses o(f) — cr 0 H(t) 
at y = ±h planes. This problem was solved by Chen [92] using finite differences as a plane- 
strain 2-D problem and has been frequently used as a reference to validate other methods. 
In our chosen geometry, height of the specimen is half the height in Chen’s specimen. This 
is done to reduce hard-disk storage requirement for the computational run on the refined 
mesh. Influence of three-dimensional effects on DSIF along the erack front is studied in 
the present work. 

Due to the symmetry in the specimen geometry and loading, only one-eighth of the 
specimen is modelled in the boundary element analysis. Fig. 4.3 shows the boundary- 
element idealization of one-eighth of the center-crack specimen corresponding to mesh no. 
4 (Table 4.1). Eight-node quadrilateral elements are used. A combination of continu- 
ous and discontinuous elements is used to reduce the overall degrees of freedom. The 
eight-node elements bordering the crack on the side of the crack surface are quarter-point 
elements, while the ones on the opposite side, i.e., ahead of the crack front are quarter- 
point, traction-singular elements. Triangular elements (collapsed eight-noded elements) 
are used to surround the intersection point of the crack front with the lateral facets. These 
triangular elements are modelled as quarter-point, traction-singular elements. 

Appropriate boundary conditions are applied to the symmetry planes. This is done by 
restraining nodal displacements in the direction normal to the plane or planes of symmetry 
to w ich the node belongs. The nodes on the crack surface, except for those directly on 
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the crack front or on the edges with lateral facets, are left free to move in three orthogonal 
directions. 

Table 4.1 shows details of some of the meshes with varying degrees of refinement. Runs 
were performed for convergence study of transient response for a straight through center 
crack under Mode-I step loading (a(t) = o Q H(t)). The profiles of dynamic stress intensity 
factor Ki(t) are seen to converge across the thickness except for the free boundary layer 
region (Figs. 4.4, 4.5, 4.6). Fig. 4.6 shows variation of DSIF across the thickness for the 
time instant tc\fh = 3.177. It is seen that mesh convergence in DSIF values is observed 
all along the crack front except for the crack vertex region. DSIFs are calculated from 
values of traction at the crack-tip nodes. The 3-D characteristics of the straight crack 
are studied using the mesh no. 4 (Table 4.1), which has 1994 nodes, 620 elements and 
5982 degrees of freedom. The selected mesh has 12 layers of elements on the face having 
the crack. The thicknesses of different layers after normalization with half-crack length 
are .008, .008, .016, .024, .032, .048, .104, .12, .12, .12, .12 and .12 respectively. Various three- 
dimensional aspects of fracture characteristics of a through thickness center crack in a 
plate under step loading are enumerated as follows: 

1. Weaker singularity at the crack vertex for Mode-I loading: For the position 
of crack vertex on the crack front, variation of normalized dynamic stress intensity 
factor (Kiit)/ (poy/im)) versus time ( tc\fh ) is shown in Fig. 4.4. At the crack vertex, 
the stress singularity is known to be weaker than square-root. Hence stress intensity 
factor, defined as the strength of the square-root singularity at the crack tip, loses its 
significance at the corner point where crack front intersects the free surface. This is 
very well reflected by lack of convergence in value of DSIF at the vertex with increas- 
ing mesh refinement. At the free surface, DSIF is much lower for mesh no. 1 and 3 
compared to mesh no. 2 as element layer adjoining the free surface is thinner in mesh 
no. 1 and 3. Value of DSIF at the crack vertex is sensitive to thickness of elements 
near the free surface. Thinner elements are needed close to free surface to estimate 
free surface effect accurately. Despite being mesh dependent, DSIF calculated from 
crack-tip tractions reflects the level of singularity at the crack vertex. 

Value of DSIF at mid-surface are close to plane-strain values as obtained by Chen [92]. 
Peak value of 2.637 occurs at the mid-surface for the specimen having b/w — 0.5. 
Value of normalized SIF for quasi-static loading at the mid-surface for the same 
specimen is evaluated as 1.088 by boundary element analysis. Parton and Boriskovsky 
[33] mentioned a 2.45 fold increase in the DSIF value over the static value in the Mode- 
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I problem considered by Chen. A factor of 2.42 is observed in the present results for 
the considered geometry. 

2. Infl uence of von Schmidt wave produced by incoming longitudinal wave 
at the free surface near the crack vertex: Variations and several local peaks 
in DSIF versus time curve caused by interaction of crack tip with Rayleigh wave, 
von Schmidt wave and reflected scattered waves for a two-dimensional plane strain 
Mode-I problem with step loading are well explained [03,94]. The present three- 
dimensional analysis brings out the influence of waves travelling back and forth in 
the thickness direction. The same cannot be observed in a two-dimensional analysis. 

A von Schmidt wave is generated by the grazing incident longitudinal wave at the 
free surface near crack vertex. This wave propagates back and forth in the thickness 
direction influencing the value of DSIF along the crack front.. Cloudier and Bishop 
[95] and Roesler [96] showed the existence of a ‘head wave' or von Schmidt wave 
trailing behind a P-wave, arising from reflection of a grazing-incident. P-wave at 
a free boundary. Beinert [97] experimentally observed the split, ting of a primary 
dilatational pulse into a train of equally spaced transverse pulses, when an elastic 
plate is loaded at its edge by a short-duration shock. Each of those transverse waves 
travels across the plate, with the velocity, c 2 , at an angle 7 given by sin 7 = c 2 /ci. 

Influence of the von Schmidt wave propagating back and forth in the thickness di- 
rection is noticeable by the mild fluctuations at regular time intervals of transit time 
a shear wave takes to travel the thickness of the specimen (Fig. 4 . 7 ). Fig. 4.8 shows 
that peak value of DSIF occurs at different positions along the crack front, at different 
instants of time. At certain instants of time, peak value of DSIF is at. mid-thickness 
position and other times peak DSIF position is shifted inside. This shows that the 
above described propagation of von Schmidt wave in thickness direction is causing 
mild variations in DSIF in the thickness direction. 

Fig. 4.9 shows the variation of divergence and %-component of curl of the displace- 
ment vector at an internal point ( Xl /w = 0 . 23 ,x 2 //i = 0 . 01 , x*/b = 0.476) with 
respect to time. Timings of the rise and fall in the curve of ‘curl’ seen in the figure 
imply the existence of a transverse wave travelling back and forth in thickness di- 
rection. But the plot of divergence does not exhibit any fluctuations corresponding 
to the times of rise and fall in the ‘curl’ curve, implying that neither Rayleigh wave 
nor dilatational wave is responsible for the mild fluctuations in values of DSIF in 
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thickness direction seen in Fig. 4.8. Only transverse waves travelling back and forth 
in thickness direction are responsible for these fluctuations. 

3. Effect of plate thickness on SIF: Influence of the von Schmidt wave is further 
highlighted in Fig. 4.10 where mid-surface Kj(t) is plotted for specimens of different 
thicknesses. Peak value of DSIF in this plot is 2.64 for b/w = 0.5, 2.41 for b/w = 0.22 
and 2.18 for b/w = 0.022. So reduction in peak value of DSIF is observed on thinning 
of the specimen due to smearing effect of greater number of traverses of shear waves 
across the thickness of the specimen. 

4. Stronger singularity at the crack vertex for Pure Mode-II loading: Plots 
of Kn(t) are shown in Fig. 4.11 for different positions along the crack front for the 
problem of symmetrically applied shear step loading on the center crack specimen. It 
is seen that DSIF is higher at the crack vertex compared to any other position on the 
crack front. Benthem [82] had found by finite difference method that a stronger than 
square-root singularity exists at the point where crack front intersects the free surface 
for a normal intersection for Mode-II elastostatic loading. Peak value of normalized 
DSIF at mid-surface position is observed to be 3.38. Influence of stronger singularity 
is reflected in the higher value of DSIF obtained at the crack vertex. Thus the present 
analysis is able to assess both weaker and stronger singularity by DSIF value at the 
crack vertex being lower or higher than corresponding plane-strain values. 

5. Effects of crack tip tunnelling: Through-thickness fatigue cracks in plate ge- 
ometries almost always show crack front tunnelling, i.e., a curved crack front with 
the point of maximum length of the crack at the center. It is also observed that 
fatigue crack fronts penetrate the free surface at a non-zero angle with the normal 
to the free surface. Bazant and Estenssoro [2] argued that this angle is fixed by the 
angle at which the square-root singularity of stresses and strains at the free surface 
penetration is restored. They showed that the front edge of a planar Mode-I crack 
that propagates in an elastic isotropic plate in a plane normal to the plate surfaces 
must terminate at the surfaces with a certain angle /3 depending only on Poisson’s 
ratio v. Tests were carried out by Rolfe and Barsom [5] on structural steel in static 
loading according to El-399 ASTM Standard. Their test results showed that the- 
oretical lines of inclination /? = 101° with the plate surface agreed very well with- 
the observed terminal directions of the first arrest marks, which corresponded to an 
essentially elastic fracture stage. Gudmundson and Ostulund [3] found that a weaker 
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singularity exists at the corner point for a dynamically growing crack and estimated 
the critical angle to be 101° almost independent of crack velocity. These studies were 
conducted using a local analysis for quarter-plane crack in a semi-infinite half-space 
subjected to symmetric loading. 

In the present study, calculations on tunnelling crack fronts are performed using cir- 
cular crack fronts with different penetration angles, viz., 95°, 98°, 101°, 105° and 
110°. The model selected for studying crack-tunnelling effect has 9 layers of elements 
of unequal thicknesses on face containing the crack. The model has 825 nodes and 
262 elements and 2475 degrees of freedom. The thicknesses of different, layers af- 
ter normalization with half-crack length are .004, .004, .008, .0 1 6, .032, .064, .128, .256 
and .358 respectively. The specimen has geometrical pa, rameters: a/w ~ 0.24, h/w = 
l,b/w = 0.22. Also in the present study, mesh is not orthogonal at the crack front, 
i.e., crack tip elements are not oriented normally to (lie crack front . The self-similar 
crack growth direction is taken as parallel to the free surface, rather t han normal to 
the crack front. Burton et al. [4] also advanced crack front, parallel to the free surface 
to produce a crack that retains its shape during propagation. 

Fig. 4.12 shows a comparison of the calculated crack-vertex A'/ (/^distributions for 
different penetration angles. It is seen that DSIF at crack vertex increases with 
increase in intersection angle. Values of Ki(t) are calculated from crack tip tractions. 
The modelled crack front is not smooth but is made of small straight line segments. 
For correct interpolation of traction singularity in crack Up elements in BEM, a 
restriction is imposed on the positions of nodes in the elements requiring that element 
edge on the crack front be straight [81]. Fig. 4.13 shows the plots of K,{1) for different 
crack front intersection angles at the mid-surface. It is seen that free surface DSIF 
is greater than mid-surface DSIF for intersection angles of 10 r and higher. At mid- 
surface, DSIF decreases with increase in intersection angle. For intersection angle of 
iOr, Fig. 4.14 shows that DSIF at crack vertex is higher than DSIF of inner layers 
implying that crack front has become critical at the corner point where crack front 
meets the free surface. 

Runs are also made for curved crack fronts on a mesh as refined as mesh no. 4 
(Table 4.1). Results from refined mesh (Fig. 4.15) show that free surface DSIF values 
are lower than mid-surface values for intersection angle of 98° and vice-versa is true 
for intersection angle of 101°. So value of critical intersection angle estimated by the 
boundary element analysis exhibits mesh-independence. Thus critical intersection 
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angle can be assessed accurately for finite domains with time-domain BEM. 

4.5 Conclusions 

1. Time-domain 3-D BEM is implemented for studying influence of free-surface on DSIF 
in a center-crack specimen under step loading. This work shows the efficacy of BEM 
and DSIF calculated from nodal traction values in studying three-dimensional aspects 
in crack/surface intersection problems for impact loading. 

2. A new partitioning scheme is proposed and implemented for spatial integration for 
time-domain BEM. 

3. Lower values of DSIF are found at corner point for Mode-I loading, confirming the 
existence of lower singularity in stress at corner point. The value of DSIF at the 
crack vertex, computed from crack-tip traction, is seen to be mesh dependent. But 
computed DSIF value reflects the influence of the exponent of singularity existing at 
the crack vertex. A lower value of DSIF is computed for a weaker singularity at the 
vertex and a higher value of DSIF is computed for a stronger singularity at the vertex. 
For the same level of mesh refinement, value of DSIF is seen to be sensitive to change 
in singularity of stresses at the vertex brought about by changes in intersection angle 
of the crack front with the free surface. 

4. Von Schmidt waves, travelling back and forth in thickness direction, influence the 
variation of DSIF along the crack front. Peak DSIF values are seen to reduce with 
the reduction in thickness of the specimen. 

5. It is possible to estimate the intersection angle at which the DSIF at the crack vertex 
is equal to or higher than DSIF of inner layers. So the intersection angle at which 
free surface vertex point becomes critical along with other parts of the crack front is 
estimated. Values of Kj(t ) at crack vertex are seen to be higher than values of K : (t) 
at inner layers for crack intersection angles of 101° and above, implying that crack 
has become critical at the vertex for intersection angle of 101°. 
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Table 4.1: Mesh characteristics for Mode- 1 runs with straight 

, crack front. 

Mesh no. 

b/w 

crack tip 
element 
size 
(M/a) 

free surface 
element 
size 
(Az/a) 

size of regular element, s 
(size/clear distance to 
crack-tip) 

(on face having crack) 

No. of 
elements 

time step 
ciAT/h 

1 

.5 

0.1 

0.004 

2.0 

340 

0.12708 

2 

.5 

0.1 

0.1 

0.5 

408 

0.06354 

3 

.5 

0.1 

0.008 

0.5 

449 

0.06354 

4 

.21 

0.1 

0.008 

0.5 

020 

0.06354 
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Figure 4.2: Center-crack tension specimen 
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Figure 4.3: Boundary Element Mesh 
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Figure 4.8: DSIF variation along the crack front at different instants of time, starting from the 
instant longitudinal wave strikes the crack face 
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Figure 4.9: Plot of xi-component of curl of displacement vector at an internal point with respect 
to time from the instant longitudinal wave strikes the crack face 













Figure 4.12: Plots of Kj(t) at crack vertex for different intersection angles of crack front with 
free surface 
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Figure 4.13: Plots of Kj{t) at mid-surface for different intersection angles of crack front with 
free surface 









Chapter 5 

Dynamic Analysis: Running Crack 


This chapter deals with the implementation of a three-dimensional time-domain boundary- 
integral formulation for a moving center-crack in a finite solid under suddenly applied 
crack face pressure. The BEM displacement time domain formulations have, hitherto, been 
limited to analyzing two-dimensional crack problems, though hypersingular formulations 
have been used to analyze finite cracks in infinite domains. In this chapter, variation 
of dynamic stress intensity factor (DSIF) along the crack front for a moving, through- 
thickness straight crack is studied for a finite solid under step loading. The state of stress 
is evaluated at the crack vertex, where crack front meets the free surface. The effect of free 
surface on DSIF is investigated. It is possible to estimate accurately the critical intersection 
angle of the crack front with the free surface at which square-root singularity is restored 
at the crack vertex under step loading for various speeds of crack propagation. It is also 
predicted in this work that for high speeds of crack propagation, crack front maintains its 
straight profile. 

5.1 Introduction 

In the present work, the standard 3-D time-domain boundary integral formulation is im- 
plemented for a self-similar moving center-crack for a finite solid under suddenly applied 
crack face pressure loading. Variation of DSIF is studied along the crack front for the 
through-thickness crack. 

A crack front intersecting the free surface in a 3-D body has been a subject of interest 
for last three decades. Through-thickness fatigue cracks in plate geometries almost always 
show crack front tunnelling, i.e., a curved crack front with the deepest point of the crack 
at the center. Benthem [1] analyzed a problem of a quarter-plane crack in a semi-infinite 
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half-space subjected to symmetric quasi-static loading, lmr a crack front normal to the free 
surface, he derived that at the corner point, stresses have r A '"' singularity for Poisson’s 
ratio (u) of 0.3. By a local finite element analysis, Ba/.ant and hstenssoro |2] estimated the 
critical intersection angle to be 101° (for v = 0.3) for a. stationary crack under elastostatic 


loading. By a finite element formulation similar to that used by Bazaut and Hstenssoro [2], 
Gudmundson and Ostulund [3] found for a half-space problem that a weaker singularity 
also exists at the corner point for a dynamically growing crack and estimated the critical 
angle to be 101°, almost independent of the crack-tip velocity. They found that, for angles 
between 90° and about 101°, singularity weakened further wit h increasing velocity and for 
angles larger than about 101°, singularity became stronger with increasing velocity. 


The success enjoyed by these local numerical analyst's is not found to the same extent in 
finite element treatments of more global or complete problems which contain crack/ surface 
intersections. Different predictions were made concerning the variation in energy release 
rate (G) with some results showing an increase in G in t he free surface region, some a 
decrease and still others no change for quasi-static analysis. Burton rt til. |4] generated 
critical crack profiles using FEM with each profile having less t han 2 % numerical variation 
in G with depth. The crack profiles generated by them were 1 significant lv less curved 
than those observed in practice [5]. Using FEM for analysis, Lin and Smith [fi] found that 
slightly non-orthogonal intersection caused a great loss of path-independence in the value of 
J integral up to 47% near the free surface. Thus in global situation, finite element analysis 
is not able to obtain an accurate stress state adjacent to the crack tip in the region close 
to the free surface, especially for cracks which are not normal to the surface. In a review 
article, Nishioka [91] has reported from an FEM analysis that value of the dynamic J- 
integral reduces near the free surfaces of the plate for a through-thickness straight running 
crack. Agrawal and Kishore demonstrated the effectiveness of HEM in evaluating state of 
stress at crack vertex and estimating critical intersection angles accurately for quasi-static 
loading [90] and for impact loading of a stationary crack [98]. 

The present work uses direct 3-D time-domain formulation of boundary element method 
to study the effect of free surface on stress intensity factor distribution along the crack front 
for a moving through-thickness crack in a finite solid under step loading. The algorithm 
used in the present work is similar to the one used by Gallego and Dominguez [10] for 
analysis of running crack in 2-D analysis. SIFs are computed along the. crack front from 
nodal values of traction in boundary element analysis [60]. Straight as well as curved 
crack profiles are studied in the present work. Critical intersection angle of crack front 
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with free surface computed in this work for step loading matches closely with the value 
of critical angle predicted by Bazant and Estenssoro [2] for elastostatic loading for a half- 
space problem. But author is not aware of any work in which critical intersection angles 
have been estimated for impact loading for a moving crack for a finite solid. 

Also for high speeds of crack propagation, a straight crack front is predicted 
by the present work. At crack vertex, DSIF value is calculated to be higher than that at 
rest of the crack front for high crack speeds, just an opposite behavior to low speed crack 
propagation or stationary crack case. Ravi-Chandar and Knauss [38] had observed with 
high speed micrography that under conditions of low stress intensity the crack front exhibits 
“thumbnail” curvature similar to the crack profiles associated with quasi-static fracture 
processes. They also observed that at increasing stress intensity the crack front straightens 
out and is identifiable as a front of multiple micro-fractures. Parton [37] experimentally 
observed that if crack velocity exceeds a certain value, then the arrival of waves reflected 
from the free boundaries to the crack tip does not affect its velocity, implying that the zone 
of process occurrence has a certain inertia, i.e., it shows some resistance to the ‘attempts’ 
to change its velocity. Based on the boundary element analysis in this study, it seems 
to the author that beyond a certain crack velocity, crack behaviour undergoes a distinct 
change or inertial terms start dominating the square-root singular terms. The possibility 
that the extent of the stress intensity factor during dynamic crack growth under transient 
conditions is more limited than a steady-state analysis would indicate was first suggested 
by the analysis of Ma and Freund [40]. Under severe transient conditions, the local crack 
tip field is not accurately represented by a stress intensity factor field. At higher crack 
speeds, the higher-order terms in the near-field asymptotic expansion can have a significant 
influence on the crack tip stress field [41] even for the case of constant crack speed, following 
fracture initiation. During the early stages of crack growth, the other terms may dominate 
the square root singular term under certain conditions. During the early phase of crack 
growth, the transient nature of stress field prevents a complete stress intensity factor from 
becoming established outside the near-tip three-dimensional zone. This feature erodes the 
value of the stress intensity factor concept in characterizing the fracture resistance of the 
material. 

A long-standing issue of fundamental importance in dynamic fracture research is the 
connection between the dynamic fracture toughness and the crack tip velocity. The debate, 
for the most part, has centered around the question of uniqueness of a relationship between 
Kf and v. Kobayashi and Dally [42], Rosakis et al. [43] and Zehnder and Rosakis [44], 
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among others, provided data sets that seem to indicate that hf v relation is reasonably 
viewed as a material property. On the other hand, the insults of Koba\ .vshi and Mall 
[45] and Ravi-Chandar and Knauss [39], based on photo-elasticity and the method of 
caustics in transmission, respectively, suggest that there is no such correspondence. The 
fundamental difficulty of achieving a well-developed stress intensify factoi field may be 
playing a significant role in these differences in results. I he pioblem may also lie in the 
assumption of linear elastic material response. Both the analysis and the models that 
underlie the interpretations of these experiments are two-dimensional in natuie. In either 
case, the three-dimensional effects that are overlooked may have too much of an influence 
to be neglected in considering the process [40]. 


For a stationary crack under suddenly applied uniform crack face pressure, Ma and 
Freund [40] found that stress intensity factor field is found only for points closer to the 
crack tip than about 5-10 percent of the distance to the cylindrical shear wave front, 
implying that the extent of the K field is so small at the onset of growth. They found that 
measured and calculated stress intensity factor histories compare very well for relatively low 
crack face pressures, but there is significant disagreement beyond crack growth initiation 
for higher pressures. The time required, before the singular stress becomes a good estimate 
of the total stress, decreases as the delay time increases. Delay time is the time between 
the onset of crack growth and the time instant when pressure is suddenly applied on the 
crack face. 


Conditions of Kf dominance exist either extremely close to the crack-tip or are even- 
tually established at long times after crack initiation. However, the higher order transient 
asymptotic representation is more successful in describing the actual field even at times 
close to the event of crack initiation or at distances relatively far away from the moving 
crack tip. The transient effect is manifested through the time derivative of the dynamic 
stress intensity factor even if the crack-tip speed is constant. 

It is necessary to apply the higher order asymptotic expansion which includes the 
transient history of the crack growth to describe the near tip deformation fields. Also 
the crack growth is indeed controlled by a material related criterion. This criterion gives 
the unique relationship between the dynamic fracture toughness I\f c and the crack-tip 
speed v. The existence of such a criterion was seen validated in a simulation by Liu and 
Rosakis [46] by using the higher order transient expansion, while the lack of the uniqueness 
of a relationship between Kf c and v has been observed when the A'j -dominant assumption 
or the steady state higher order expansion is used [39,45]. 
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In the present study, the author has investigated if it is possible to capture this loss of 
square-root dominance by the boundary element analysis. SIF is calculated from crack-tip 
or crack front tractions in BEM. Value of traction at a node on crack front exhibits an 
integrated influence of level of singularity over the elements adjoining it. Any weakening 
of singularity, for example at crack vertex due to free surface effect or loss 
of dominance of square-root terms at high speeds of crack propagation, gets 
reflected in the value of SIF calculated from nodal tractions. 

In the past three decades, the application of the finite element method (FEM) to 
dynamic crack propagation has made important advances. Interesting reviews can be 
found in the papers by Aoki et al. [47] and by Atluri and Nishioka [48]. Two different 
approaches followed are based on one of the following two concepts: a stationary mesh, 
which includes a ‘node release’ mechanism, or a moving mesh which normally includes a 
singular element. 

The BEM has appeared as an alternative for elastic fracture mechanics problems. This 
method seems to be a better choice than FEM for elastodynamic fracture mechanics be- 
cause the discretization is restricted to the boundary surface and the concept of singular 
element is simplified. In particular when dealing with dynamic crack propagation, the 
remeshing is conceptually much simpler in BEM than in any domain technique. Classical 
methods like finite difference method or FEM are difficult to apply to the study of crack 
growth because of their intrinsic numerical dissipation of high frequency elastic waves. 

Several different boundary element formulations exist for study of dynamic fracture. 
Hypersingular formulations have been used to analyze finite plane cracks in 3-D infinite 
domains [83,84]. Laplace transform BEM has been used to analyze cracks in finite solids un- 
der dynamic loading [85]. Wen, Aliabadi and Rooke [86,87] presented a Laplace transform 
displacement discontinuity and fictitious stress formulations for 2-D and 3-D problems. 
They determined the stress intensity factors for several 2-D and 3-D problems using an 
equivalent stress approach. Das [99] used the displacement discontinuity method to deter- 
mine the slip of a planar crack propagating in an infinite medium. Das and Kostrov [100] 
presented displacement and traction equations for 3-D dynamic shear crack propagation 
problems. Roller et al. [101] used the traction equation expressed in terms of displace- 
ment discontinuity to analyse dynamic anti-plane cracks. Fedelinski el al. [102] extended 
the dual boundary element formulation in time-domain to analyses of mixed- mode fast, 
growing cracks. Time domain formulation is advantageous where an accurate solution 
is required from the very beginning and is the best possible approach for problems with: 
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time-dependent geometry [50]. 

Dominguez and Gallego [60] used the din'd, time domain boundary element formulation 
in combination with the traction-singular quarter -point element:, that had been used by 
Martinez and Dominguez [73] and Chirino and Dominguez jss' tor static and time 'harmonic 
computations, respectively. Gallego and Dominguez jll); modelled the ,-tack growth in 
two-dimensional analysis by using moving singular (dements and a reme.shing technique 
For moving elements, they used space- and time-dependent intei point ion functions The 
same algorithm is extended in the present work to 3-D analysis of problems of dynamic 
crack propagation in finite bodies. To author's knowledge this is the first time that direct 
displacement time-domain 3-D BEM is used to study dynamic cta.i propagation in finite 
bodies using moving mesh procedure. 


5.2 Analytical problem of .self-similar, constant, veloc- 
ity, crack-propagation from a finite initial length 

A crack propagation problem is described Imre to facilitate e, tson with rial 
resuits later. The analytical problem considered is tl„. g.mvth ,,t a half plane crack in 
an otherwise unbounded body of a homogeneous. i»„ru,„e ebs„e s„l„i Suppose that a 
uniform normal pressure of magnitude «•„ l.egins to act on I he faces „f 1 1,„ nark at a certain 
mstant, say time i = 0. The crack starts growing after an initial delav , , . The vmmt 

expands over the newly-created era* fac.es. The analytical solut,,.,, f.„ dvnatnie stress 

mtenstty factor in the interval 0 < ( < T ( U , the crack is stationary l e. 


K(t) = «n,V,Vr 


( 5 . 1 ) 


“ao * fT TOV ° S,,Ml f0r U '" ""*■» «»■ » « fmo t. I’oisson’s 

crack is ' * > r ’ UlC dy,,amir »*«* "ntv hr... -o, the propagating 
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where the slowness function S+{l/v) is an integral arising from Wiener-Hopf factorization 
for the moving crack problem by Freund [11,12]. Freund [12] notes that S + ( 1/v) is “not 
too different from unity over the full range of its argument” and it is given as 


S+(l/v) = exp[- 


r^i 


tan 


7T 


I /n n / ON 2 J 


dc 


r] (5-4) 


* C'2 


(2 - c 2 /4) 2 J c(c - v) 

Values for S+( 1/v) and k(v) for Poisson’s ratio of 0.3 for several values of v/cr are shown 
in Table 5.1. Here, B is defined as 


B = 2 vTEMd (5.5) 

1-v 

For the particular situation of crack face pressure expanding over the newly created crack 
faces, DSIF is independent of arbitrary delay time after t> r. 


5.3 Transient Boundary Integral Formulation 

Consider a homogeneous, isotropic and linearly elastic body of volume V bounded by a 
regular surface S, with its mass density p and dilatational and shear wave speeds as Ci and 
c 2 respectively. Using Stoke’s state of quiescent past corresponding to a point load with a 
Dirac delta function variation, Love’s integral representation can be derived using Betti s 
reciprocal theorem. For zero initial conditions and zero body forces, the displacement of a 
point y can be represented by the following boundary integral equation: 


Cij(y)ui{y,T ) = 

J [Gij(x, y , T) * ti{x,T ) - Hij{x,y, T) * Ui(x, T)} d S(x) 

(5.6) 

where 

rT 



Gij*ti= / Gij{x,T;y,T)ti{x,T) dr 

J 0 

(5.7) 


Hi, * Ui = 1 Hij(x,T-,y,T) Ui{x,r) dr 

Jo 

(5.8) 


are Ricmann convolution integrals. In the above, «,(*, r) and U(x, r) are the displacement 
and traction components respectively at point x and time r. The free-term tensor Cvj is 
calculated by general close-form expression developed by Mantic [56]. The fundamental 
solutions Gy and JEZ# are the displacements and tractions in direction i at a point z at 
time T due to a unit impulse vector in direction j applied at a point y at a preceding time 
r . The fundamental solutions are given in explicit form as 
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+a ii {(l/ci)5(v - r/c\) - ( l/cl)5{v - r/c 2 )} 4- r/r 2 ) 


(5.9) 


where v = T-r ; a# = Vr 3 ; % = %/r; hi = x, - ?/i ; r is the distance between 
points x and y. Also % is the Kronecker’s delta and <J is the Dirac delta function. 


H iS (x,T-,y, r) 


J_r 

47T 


/•1/C2 

- 6*1(5% - ey) / A<$(v - Ar)dA 4- (12% - 2c, ; ) {<)'(» - r/c 2 ) 

di/d 


-{c 2 /ci) 2 6(v - r/ci)} 4- 2rdij/c 2 {5(v - r/c 2 ) - {c 2 /c\YS(v - r/r { )} 
- %•( 1 - 2c|/cf){(5(t; - r/ci) + ( r/a)8(v - r/c x )} - Sy{<5(*' - r/c 2 ) 
+(r/c 2 )5(v - r/c 2 )} 


(5.10) 


where v — T t; /i, — x, t/g — hihjh m ‘ii m /i , fij •- hjHj/i , </q ( hiiij 4- 

$ijh m n m )/r 3 ; = % 4- gij . Also n* are the components of outward pointing normal 

vector and n m h m = nj/ii 4- n 2 h 2 4- n 3 /i 3 . 


5.4 Numerical implementation 

The time-domain BEM in three-dimensions consists of two basic steps as explained in 
[49,58,59]: 

1. A discretization of the real time axis into a sequence of equally spaced time intervals 
with the assumption of linear variation of displacements and tractions over each 
time interval. When piecewise linear variation is used for both displacements and 
tractions, the time-stepping scheme is supposed to become unstable [GO, 61]. 

But in this study, piecewise constant interpolation for traction variable is seen to 
give poor results in 3-D time domain analysis. So traction variable is also linearly 
interpolated in time like displacement variable. Good results are obtained and insta- 
blity is not seen to affect the results. It is, in contrast, with stationary crack results 
where linear time-stepping for traction variable results in. instability. It is possible 
that for moving elements, self-effect dominates over the other-effect [103] for lin- 
ear interpolation of tractions over time and hence, resulting in stablity of numerical 
scheme. 


2. A discretization of boundary into 8-noded elements. 
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The time span of interest is divided into N equal time steps of size AT so that 

T n = n AT, n = 1, 2, . . . , N. (5.11) 

On the basis ol these discretizations a time stepping solution of equation (5.6) can be 
established for the boundary displacements and tractions over each element and for each 
time step. 

The displacement boundary integral equation can be written as 

r'l'N r rT^-i r 

CijUi(y, T n )~ / [Gijti - H ijUi } dSdr = / / [G tJ U - H ijUi ] dSdr (5.12) 

Js Jo Js 

where the integral on the right hand side is the contribution due to past dynamic history. 
This equation is also an exact statement since no approximation has yet been introduced. 
However, in order to solve the above equation, one has to approximate the time variation 
of the field quantities, in addition to the usual approximation of spatial variation. The 
time integration in the above equation is done analytically, and the spatial integration is 
done numerically. A system of algebraic equations is obtained at time step N through 
nodal collocation, i.c., by allowing y to coincide sequentially with all the nodal points of 
the boundary. 

Equation (5.12) can be written in a more condensed form as a typical BE equation: 

= EEGtTc < 5 - 13 ) 

n=l (?=1 n=l q — 1 

or, in matrix form: 

N N 

E 

71=1 n - 1 

where u 11 and t n include the three displacement and traction components of all the 
nodes at time step n. Here, p and q denote the collocation point and integration element 
respectively. And Q are the number of elements in the boundary element mesh. 


H Nn u n 


= Eg 


NrijXi 


(5.14) 


5.4.1 Temporal and spatial interpolation of displacement and 
traction 

The boundary variables are interpolated in terms of their values at nodes by 

«j(x,r) = y j [Mj + MpUi(x)] 

71=1 


(5.15) 
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ti(x,r) = y; [M] 1 t*~ l (x) + i\/,” /"(a:)] (5.16) 

71=1 

and 

«?(*) = [{(*)]«?', <-1.2,3 (5-17) 

C=1 

urn*) = i==I > 2 > 3 ( 5 - lg ) 

C=1 

where uf c and «| n_1)c represent nodal values of displacement at time nodes n and (rc-1) 
respectively. Similar interpolation is done for traction variables also. Here, N is the total 
number of time steps. Mr and Mp are the temporal interpolation functions, related to 
local time nodes I (initial) and F (final), and are of the following form for piece-wise linear 
variation: 

U}‘ = - Mr) (5.19) 

Mr = ( 5 ' 2 °) 

where 

, f \ _ j 1 if (n - 1)A T < r < 7i AT 

n T \ 0 otherwise 

or 

<pn{r) = H (r - (n - 1)AT) - H{r - nAT) (5.21) 

H being the Heaviside function. 

Both displacements and tractions are interpolated by piecewise linear functions in the 
present work. 

For spatial interpolation of boundary variables, an approach [62] which combines the 
advantages of continuous and discontinuous elements in a single formulation is used in 
this work. In this formulation, expressions of the shape and interpolation functions are 
unified in a form suitable for fully continuous, fully discontinuous or transition elements 
(i.e., discontinuous only along one or more edges). Three types of boundary nodes are 
identified depending on the function they fulfil: 

1. Geometrical nodes: those used for interpolating the boundary geometry. 

2. Interpolation nodes: the points on the boundary which define the interpolation of 
the variables. The system nodal values correspond to these nodes. 
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3. Collocation nodes: where the boundary integral equation is enforced. 

One property of the boundary element formulations is that the discontinuities of the vari- 
able across the element edges do not invalidate the convergence of the technique [63]. 
Collocation and interpolation nodes for continuous elements coincide with the geometrical 
nodes. Collocation nodes and interpolation nodes for discontinuous elements and discontin- 
uous side of transition elements are situated inside the element. For quarter-point element, 
geometrical and interpolation nodes coincide, but collocation corresponding to nodes on 
crack front is done at nodes situated inside the elements on either side of the crack front. 
To reduce the number of equations and enhance accuracy, equations corresponding to 
collocation nodes before and after a given interpolation node are added up [64]. 


5.4.2 Analytical integration of temporal terms 


The time integration in equation (5.12) after using equations (5.15) and (5.16) is done 
analytically and the spatial integration is performed numerically. Explicit expressions for 
time integrals may be seen in [49,58]. Temporal integration boils down to integration of 
terms like / Q Ay 5 (T — t — r/c) Mj (r) dr, J AT X5 (T — r — A r) Mj (r) dAdr etc. 

For example, 




(1— r/AT)0i(r) (r) dr 


(l-^) + Sf r if c(T - AT) < r < cT 
0 otherwise 

(5.22) 


Similarly, all other temporal terms can be expressed in the form a + br + £ after analytical 
integration in time. For moving elements, temporal integration also results in an additional 
r 2 term. It is noticed that for the distance reached by both shear and compression waves 
travelling from a source node up to a time AT, coefficient of A? is zero. Thus, the temporal 
terms, i.e., time integrations are not adding to the singularity in r. Here, c is compression 
(ci) or shear ( c 2 ) wave speed. 


5.4.3 Spatial integration 

The kernels in elastodynamic case have discontinuous derivatives at points such that 
r/ci( 2 ) = kAT {k integer). Therefore, the standard numerical methods, like Gaussian 
integration, give poor results if applied over the whole element. So a partitioning scheme 
is used for elements so as to facilitate accurate spatial integration [98]. 

Principal value singularity occurs for traction kernel Hij when source element becomes 
also a receiver. The Cauchy principal value integral is treated in the same way, as done 
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by Guiggiani and Gigante [7] for the elastostatic. case with the help of partitioning scheme 
proposed in [98]. The non-linear transformation suggested by Doblare and Gracia [66] is 
also incorporated. 

Angular variable transformation of Hayami [79] is used to weaken the angular near- 
singularity arising when collocation node is near the edge of an element.. 

5.4.4 Crack Elements 

Quarter-point elements are used along the crack front and traction-singularity is incorpo- 
rated in the quarter-point elements ahead of the crack front [71,72]. The stress intensity 
factors can be computed directly by the boundary element code as follows: 

I(j = 

Kjj = t\'/2'nl (5.23) 

Km = il^Txl 

where U are the tractions at nodes on the crack front and l is the length of crack front 
elements in the direction of singularity. It may be noted that, the use of traction nodal 
values of the singular element at the crack tip for SIP computation is substantially less 
sensitive to the discretization than any of the displacement, correlation procedures [60]. 


5.5 Formulation for a moving singular element 

This formulation is useful for multi-domain or symmetric modelling to simulate crack 
growth. In multi-domain modelling, a boundary is introduced along the crack, following the 
known direction of propagation. Because of symmetry of the problem under consideration, 
only one-eighth of the domain is analysed. Assume that crack propagates at, a speed v 
and that at time T , the position of the crack front and the discretization of the part, of 
the boundary that contains the crack are shown in Fig. 5.1. The elements adjoining the 
crack front are quarter-point elements. After a time increment AT, the crack front has 
moved vAT. As the crack front advanced from left to right, all the elements before and 
after the crack front are also shifted to the right by vAT. Only centroids of the leftmost 
and the rightmost columns of elements shift by |t>AT, the leftmost and the rightmost 
edges remaining stationary. Thus, the leftmost column of elements increase in size and the 
rightmost column of elements decrease in size and elements in between these two columns 
of elements remain the same size. Crack propagation is studied in the present work till the 
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time the rightmost column of elements can be decreased no further in size. To start with a 
mesh is prepared having small elements in the leftmost column and large elements in the 
rightmost column, so as to allow large amount of crack growth without requiring addition 
of new column of elements to the left of crack front and deletion of an existing column of 
elements to the right of crack front. Above scheme of element movement helps in reducing 
computations of matrices, as will be explained later. In the boundary element mesh shown 
in Fig. 4.3, elements on the bottom plane move with the crack front. Elements on the 
front, top and the back faces are also moved with the crack front in a similar fashion as 
the elements on the bottom face. 

The BE matrix for time step N can now be written in a subdivided form as 

, „ x N 


N 

E 

n=l 


Nn 

fi 


H 

H Nn 
of 


rrNn 

H fc 

rrNn 
£1 cc 


= E 

n= 1 


° 7 / 

G Nn 
cf 


G Nn 

fc 

G Nn 

rr 


(5.24) 


where the subindex ‘f ’ stands for the nodes whose position remains fixed during the remesh- 
ing process and ‘c’ for those that change position during that process. It is important that 
when equation (5.24) is written for a new time step N, the sub-matrices corresponding to 
changing collocation nodes and/or to changing collocation elements have to be calculated 
anew for all the previous time steps n. The translation property of the fundamental so- 
lution cannot be used in those cases because the position of the nodes at time step N is 
not the same as before. Thus H* npq and G% npq depend on the N and n values and not 
only on the difference N — n. But for the combination of moving elements and moving 
collocation nodes, where there is no change in relative position of collocation nodes and 
integration elements, translation property can still be made use of, starting from the time 
interval when crack starts propagating. All moving elements, except for the leftmost and 
the rightmost columns of elements, fall into this category. 

For moving elements, the space interpolation functions of the variables move with 
the elements and therefore, these functions are not only space dependent but also time 
dependent. Displacements and tractions for moving elements are represented as follows: 


Ui(x,r) = ^[M r n uf(x,r) + M$u?{x,t)] (5.25) 

n= 1 

U(x,r) = J2[M?t!{x,r) + M£t({x,T)] (5-26) 

n= 1 

la the regular formulation of the time domain BEM, the space shape functions are separated 
from the time integrals which are done analytically. This is not possible when space shape 
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functions are time-dependent. In order to do the time intcgiation, an appioximation of the 
space shape functions’ time dependence is used [10], Two terms of their series expansion 
are taken: 

= u?(x) + (t -T„)v"(x) 
t!(x, r) = ^- 1 ( a! ) + (T-T„.,)«r 1 W 

= «?(*) + (t-T„)(?(x) (5.27) 

Assume that for each time step the elements move with a constant, speed v in x , -direction 
same as that of the crack front, the shape functions can be written as 

Ui(Xi,X 2 ,X 3 ,T) = Ui(X{ - VT,X2 ,Xs) 

= Ui(x,x 2 ,x: 0 (5.28) 


The functions on the opposite sides of the equation in (5.28) are diifereut, of course, but 
the values of the functions represent one and the same physical quantify. The ,r, x 2 , x 3 - 
coordinate system translates at speed v in the aq- or x-direotion with the origin of coordi- 
nates in this system on the crack front. 

Then, 


u i (xi,X2,x 3 ,r)| T=Tn _ 1 


dvfa 1 (x) 
V fa 


(5.29) 


Ui(xi,X2,X3,T)\ r -T n 


du r Mx) 

V — ~ — 


(5.30) 


The expressions of shape functions for iq(x) and ti(x) are those of the usual shape functions 
and their derivatives are easily obtained. In the present work for traction singular elements, 
ti(x) is singular of the kind x _1/2 which makes its derivative singular of the kind 
and the concerned integral a hypersingular one. Gallego and Dominguez [10] computed 
the finite part of these integrals in their 2-D analysis of moving crack problem. Similar 
approach is followed and finite part of these integrals is computed as explained later. 


5.6 Finite-part integration 

Owing to the unknown transformation properties of the finite-part integrals undergoing 
a non-linear coordinate transformation, the definition formula of the finite-part integrals 
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is applied in extrinsic coordinate system. The resulting integrals are regular an d may be 
evaluated by standard Gaussian quadrature rules [104]. 

In this work, integration is performed in extrinsic polar coordinates with origin at 
projection of collocation node on the plane tangent on the closest point on the element 
under integration (Fig. 5.3). Corresponding to a Gauss point G , a point G a is found on 
the crackfront along the direction of singularity of stresses. Subtraction of value of the 
kernel at point G 0 from the kernel at point G makes the integrand regular by weakening 
the x~ z ! 2 singularity. Corresponding term which is added requires finite part integration. 
Corresponding to (p, 9) coordinates of point G , intrinsic coordinates (£, rf) are calculated, 
thus, facilitating evaluation of shape functions at the Gauss point G. 

Finite-part integration over the element is done by choosing Gauss points along the 
crack front and doing finite-part analytical integration along the direction of singularity of 
stresses for each Gauss point. 


5.7 Details of moving mesh algorithm 


Equation (5.16) for the interpolation of traction variable U for eight-noded elements can 
be written as 


N 8 


U{x,r) = 


n=l c= 1 


iV c [£(:c)] — v(r — T n _i) 


dN c [£(x)} 


+ 


N c [Z(x)]-v(T-T n ) 


dN^x)} 


dx 


dx 


r=T n 


M J (r)4 n_1)c 


r~T n ~i 

M e (r)C 


(5.31) 


Similar expression can be written for the displacement variable. In the chosen problem, 
crack front is advancing in X\- or ^-direction. Hence, derivative of shape functions at 
instant t = T n _i and r = T n are calculated from the geometry of the element at time 
nodes T„_i and T n respectively as given below: 


dN c 


dx 


r-Tn - 1 Or r=Tn 


dN c d£ 
d£ dx 


+ 


r=T„_ i or r=T n 


dW dv 

dr] dx 


( 5 . 32 ) 


T—T n ~\ or r=r n 


By substitution of Equation (5.31) into the boundary 
integrals can be separated from the space integrals and 
be written as 
r rTn 


/ / Gij{x,T^; y,r)ti(x,r)drdS 


integral equation (5.6), the time 
the coefficients of the system can 
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rr»T n 


In above integral over element 5 P , temporal integrations are done analytically. Dis- 
placement interpolation in a quarter-point element with crack front at. £ = - 1 (say) can 
be expressed as 


«i(f,»7) = a*(>7)f 2 + &<(»?)£ + <*(*/) 

= ai{r])(2y/xjL - if + hi\‘i])(y/x/L - 1) 1 <',,(//) 

= ^(4oi(»7)) + ^-(-4 Oi(r/) + &<(»/)) + («,.(//) - 6i(r/) -1- r f (»/)) (5.34) 

where terms s and L denote the distance of field point and extent of crack element, respec- 
tively from the crack front along the direction of singularity as explained for quadrilateral 
elements in Fig. 5.3. Traction interpolation in a traction-singular quarter-point, element 
with crack front at £ = — 1 can be expressed as 

Ufay) = ^ [“iivK 2 + + Ci(rj)] 

= + (— 4aj(j7) + - bi(?/) + 

~ Mv) + B i( r l)\jj i + (5-35) 

The derivative of traction variable with respect to coordinate along the direction of 
crack propagation x (considered here) can be written as 

I = Hvb- C ^2 + 

( d Mv) . dBiiri) fx dCi{r]) fl\ dr / 

\~w~ + ~dru + -drn)Fx 

\fl + C "W][^+D u ^j- 2 


= ^.ii(r?) + B u {ri) 


(5.36) 
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where a has value 0° for quadrilateral elements where quadrilateral sides having quarter- 
point nodes are parallel to the direction of crack propagation. For triangular elements at 
the crack vertex in the modelled geometry, a represents the angle made by the ray joi nin g 
the field point to the crack vertex with the direction of crack propagation. Here, the term 
containing is integrated by Gauss-Jacobi integration rules. The term containing 
requires finite-part integration. Integration is performed in extrinsic polar coordinates p 
and 0 and is treated as follows: 

[ J F (P’ 6 )^P d P de = J J( F (P(^V)^{CV))- F o{^= -l,r?))^Jpdpd0 

+ ~^f2 dxdr] ( 5 -37) 

where J(r]) is the Jacobian for mapping of the element edge along the crack front to 
intrinsic coordinates rj = — 1 to 1. Here, L and x can be written in terms of p and 6 for 
straight-edge elements in the form (a{6)p + b(Q)). First integral on the right hand side can 
now be calculated by Gauss-Jacobi integration rules. And the second integral on the right 
hand side is written as 

j' ^ J(v) F o(£, = -M) ^ -^ dxd v = Y^WiJ{Vt) F o{£ = -i,v = Vi) F - R J o ~^ dx 

And 

F - R [ L ^ = - 2 ( 5 - 38 > 

5.8 Results and Discussion 

Elastic analysis of a running crack subjected to uniform step pressure on crack faces has 
been performed for this work. Fig. 5.4 shows the geometry of the center crack specimen 
used in the numerical analysis. The specimen has the parameters: a = 0.24 cm, w = 
1 cm, h = 1 cm and b = 0.21 cm where a is the crack half-length, w is specimen half-width, 
h is half-height and b is specimen half-thickness. The material is assumed to be linear 
elastic with properties as: shear modulus p = 90 GPa, Poisson’s ratio v = 0.3 and density 
p = 7800 kg/m 3 . Hence, c x = 6354 m/s and c 2 = 3396 m/s. The center crack specimen is 
loaded by suddenly applied uniform normal pressure to crack faces cr yy U) = —cr 0 H(t). This 
problem was solved analytically by Freund [11] for a semi-infinite crack in an unbounded 
domain. Analytical solution [11] is valid until unloading waves from the free boundary or 
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waves diffracted from the other crack tip reach the crack tip under consideration. Influence 
of three-dimensional effects on DSIF along the crack front is studied in the present work. 

Due to the symmetry in the specimen geometry and loading, only one-eighth of the 
specimen is modelled in the boundary element analysis. Fig. 4.3 shows the boundary- 
element idealization of one-eighth of the center-crack specimen. Eight-node quadrilateral 
elements are used. A combination of continuous and discontinuous elements is used to 
reduce the overall degrees of freedom. The eight-node elements bordering the crack on the 
side of the crack surface are quarter-point elements, while the ones on the opposite side, i.e., 
ahead of the crack front are quarter-point, traction-singular elements. Triangular elements 
(collapsed eight-noded elements) are used to surround the intersection point of the crack 
front with the lateral facets. These triangular elements are modelled as quarter-point, 
traction-singular elements. 

Appropriate boundary conditions are applied to the symmetry planes. This is done by 
restraining nodal displacements in the direction normal to the plane or plain's of symmetry 
to which the node belongs. The nodes on the crack surface, except, for those directly on 
the crack front or on the edges with lateral facets, are left free to move in three orthogonal 
directions. 

Table 5.2 shows details of the mesh used for boundary element analysis. The mesh 
has 1994 nodes, 620 elements and 5982 degrees of freedom. It has 12 layers of elements 
on the face having the crack. The thicknesses of different layers alter normalization with 
half-crack length are .008, .008, .016, .024, .032, .048, .104, .12, .12, .12, .12 and .12 respectively. 
Various three-dimensional aspects of fracture characteristics of a through thickness center 
crack in a plate under step loading on crack faces are enumerated as follows: 

1. Weaker singularity at the crack vertex for Mode-I loading at low crack 
speeds: Fig. 5.5 shows the DSIF variation with time for the case of a stationary 
crack under suddenly applied crack face pressure. The plot clearly shows the influence 
of lower stress singularity at crack vertex resulting in lower DSIF values at the crack 
vertex. At the crack vertex for a stationary crack, the stress singularity is known 
to be weaker than square-root. In these plots, z/b = 0.0 denotes the free surface 
position along the crack front and z/b = 1.0 denotes the mid-surface position. 

Figs. 5.6, 5.7 show the computed values of normalized Mode-I DSIF versus normal- 
ized time for crack propagation velocities v = 0.2c fi and 0.4c ft , respectively. The 
boundary element results are compared with analytical solutions for 0.2c/{ and 0.4ck, 
denoted in the figures by dotted lines. It can be seen from the figures that the 
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agreement between the boundary element mid-surface DSIF values and plane-strain 
analytical results is quite good for v = 0.2 cr. The mid-surface curve for v = OAcr 
lies quite below the plane strain analytical curve. This reflects the effect of loss of 
dominance of square-root singularity at high speeds of crack propagation. The an- 
alytical solutions are for the case of half-plane crack propagating in an unbounded 
domain under the action of suddenly applied uniform normal pressure acting on the 
crack faces as discussed in section 5.2. 

For a running crack, reduction in value of DSIF at crack vertex is noticed at low 
crack speeds (Figs. 5. 6, 5. 7). At crack vertex, value of DSIF is lower than the rest of 
the crack front for v/cr — 0.2 and 0.4. But this behaviour changes at speeds beyond 
v/cr = 0.4 (Figs. 5. 8, 5.9). For higher crack speeds, DSIF value at crack vertex 
becomes higher than the rest of crack front. This may be the result of greater erosion 
of dominance of square-root singularity along the crack front inside the specimen 
rather than near the crack vertex at higher speeds of crack propagation. Also this 
implies that free surface retardation of crack front would not occur at high crack 
speeds. This result is in tune with experimental observations where at high crack 
speeds, straight crack front is observed. Ravi-Chandar and Knauss [38] explained 
that there is no longer a single crack front, but an ensemble of cracks propagating 
along the apparent crack front. 

Figs. 5.10,5.11 show plots of DSIF along the crack front at various instants of time for 
two different crack velocities respectively. It is seen that the curves have oscillations 
in their values along the crack front. Mid-side nocles have lower values of DSIF 
than corner nodes of elements on crack front. This may be a result of sensitivity of 
boundary element method to spatial instability [101], which is manifesting itself in 
crack propagation problems. Figs. 5.12,5.13 show DSIF variation along crack front 
by considering tractions at either corner or mid-side nodes for DSIF computations 
respectively. These plots clearly show increased DSIF value at crack vertex for crack 
speed of v/cr = 0.8. 

Again, Figs. 5.14,5.15 show the variation of DSIF with time at the crack vertex and 
mid-surface respectively for time step c x AT/h = 0.06354 for different crack speeds. 
Figs. 5.16,5.17 show the same plots for a reduced time step size ciAT/h = 0.015885. 
It is clearly seen that around or beyond v/c R = 0.4, free surface values of DSIF are 
higher than mid-surface values, implying that there will be no crack retardation at 
free surface at speeds beyond v/c R = 0.4. The results seem to indicate that that be- 
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yond a certain crack velocity, crack behaviour undergoes a distinct, change 01 inertial 
terms start dominating the square-root singular terms. Parton [37] experimentally 
observed that if crack velocity exceeds a certain value, then the ai rival ol waves re- 
flected from the free boundaries to the crack tip doers not affect its velocity, implying 
that the zone of process occurrence has a certain inertia, i.e., it shows some resistance 
to the ‘attempts’ to change its velocity. So crack speed v = 0.4c/; appears to be the 
‘critical’ speed beyond which behaviour at crack vertex gets reversed. The value of 
DSIF at crack vertex becomes higher than the rest of the crack front, just opposite 
to the behaviour for a stationary or low speed crack. 

Another observation, worth noting, from these figures is that stationary part of DSIF 
curves are much below the stationary analytical solution for smaller time step plots. 
This concurs with the conclusion of Ma and Freund [40] that time is needed for the 
stress intensity factor controlled field (A'f-dominant field) to be fully established. 
This shows that DSIF values computed from boundary element, analysis are sensitive 
to loss of square-root dominance in transient analysis. 

2. Effects of crack tip tunnelling: Through-thickness fatigue cracks in plate ge- 
ometries almost always show crack front tunnelling, i.e., a curved crack front, with 
the the point of maximum length of the crack at the center. It, is also observed that 
fatigue crack fronts penetrate the free surface at a non-zero angle with the normal 
to the free surface. Bazant and Estenssoro [2] argued that this angle is fixed by the 
angle at which the square-root singularity of stresses and strains at the free surface 
penetration is restored. They showed that the front edge of a planar Mode-I crack 
that propagates in an elastic isotropic plate in a plane normal to the plate surfaces 
must terminate at the surfaces with a certain angle ft (Fig. 5.18) depending only on 
Poisson’s ratio v. Tests were carried out by Rolfe and Barsom [5] on structural steel 
in static loading according to E-399 ASTM Standard. Their test results showed that 
theoretical lines of inclination ft = 101° with the plate surface agreed very well with 
the observed terminal directions of the first arrest marks, which corresponded to an 
essentially elastic fracture stage. Gudmundson and Ostulund [3] found that a weaker 
singularity exists at the corner point for a dynamically growing crack and estimated 
the critical angle to be 101° almost independent of crack velocity. These studies were 
conducted using a local analysis for quarter-plane crack in a semi-infinite half-space 
subjected to symmetric loading. The present study deals with the analysis of a center 
crack finite solid under suddenly applied uniform normal pressure to crack faces. 
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In the present study, calculations on tunnelling crack fronts are performed using 
circular crack fronts with different penetration angles, viz., 95°, 98°, 101°, 105° and 
110° using a mesh similar to that of Table 5.2. Also in the present study, mesh is 
not orthogonal at the crack front, i.e., crack tip elements are not oriented normally 
to the crack front. The self-similar crack growth direction is taken as parallel to the 
free surface, rather than normal to the crack front (Fig. 5.18). Burton et al. [4] also 
advanced crack front parallel to the free surface to produce a crack that retains its 
shape during propagation. 

Fig. 5.19 shows a comparison of the calculated crack- vertex Kj{t) distributions for 
different penetration angles for an example case of crack speed v = 0.4c^. It is 
seen that DSIF at crack vertex increases with increase in intersection angle. Values 
of Kj (/.) are calculated from crack tip tractions. The modelled crack front is not 
smooth but is made of small straight line segments. For correct interpolation of 
traction singularity in crack tip elements in BEM, a restriction is imposed on the 
positions of nodes in the elements requiring that element edge on the crack front be 
straight [81]. Fig. 5.20 shows the plots of Ki{t ) for different crack front intersection 
angles at the mid-surface position. It is seen that free surface DSIF is greater than 
mid-surface DSIF for intersection angles of 101° and higher. At mid-surface, DSIF 
decreases with. increase in intersection angle. For intersection angle of 101°, Fig. 5.21 
shows that DSIF at crack vertex is higher than DSIF of inner layers implying that 
crack front has become critical at the corner point where crack front meets the free 
surface. Thus critical intersection angle can be assessed accurately for finite domains 
with time-domain BEM. 


5.9 Conclusions 

1. Time-domain 3-D BEM is implemented for studying influence of free-surface on DSIF 
in a center-crack specimen under suddenly applied uniform pressure to crack faces. 
This work shows the efficacy of BEM and DSIF calculated from nodal traction val- 
ues in studying three-dimensional aspects in crack/surface intersection problems for 
impact loading. 

2. A moving mesh procedure is proposed and implemented for 3-D time-domain dis- 
placement boundary element analysis. 
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3. Lower values of DSIF are found at corner point for Mode-I loading for low speeds of 
crack propagation, confirming the existence of lower singularity in stresses at corner 
point. But at crack speeds beyond v = 0.4cr, higher values of DSIF are computed 
compared to the rest of the crack front. The study correctly predicts that for high 
speeds of crack propagation, crack front maintains its straight profile. 

4. Values of DSIF calculated from nodal tractions at crack front show sensitivity to loss 
of square-root dominance in transient fields. 

5. It is possible to estimate the intersection angle at which the DSIF at the crack vertex 
is equal to or higher than DSIF of inner layers. So the intersection angle at which 
free surface vertex point becomes critical along with other parts of the crack front is 
estimated. Values of Ki(t) at crack vertex are seen to be higher than values of K,(t) 
at inner layers for crack intersection angles of 101° and above, implying that, crack 
has become critical at the vertex for intersection angles of 101". 
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Table 5.1: Values of the Freund [11,12] slowness function S+(l/v) for various crack speeds v for 
Poisson ratio v = 0.3 


V/C R 

S + (l/v) 

k(v) 

0.2 

0.96032 

0.87766 

0.4 

0.90746 

0.73836 

0.5 

0.90746 

0.65976 

0.6 

0.90746 

0.57280 

0.8 

0.71753 

0.35859 
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Table 5.2: Mesh characteristics fo r Mode-1 runs with strair.ht 


Mesh no. b/w 

crack tip 

free surface 

sizo oi regular elements 


element 

element 

(size/ dear distance to 


size 

size 

crack-tip) 


{M/a) 

( Az/a) 

(on lace having crack) 

1 .21 

0.1 

0.008 

0.5 


nark trout.. 

of 1 imc step 
cli'tin ’tit s I’lST/h 


ti’.'tl 


0.00354 



5.9 Conclusions 


15] 






152 


Dynamic An alysis: Running Crack 



Figure 5.2: Boundary Element Mesh 





Figure 5.3: Quarter-point element projected in the plane tangent to the nearest node. 
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Figure 5.5: Plots of Kj{t) at different positions on crack front for suddenly applied crack 
face pressure and comparison with analytical solution of problem of semi-infinite crack m an 
unbounded body 
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Figure 5.6: DSIF variation with time for different positions along the crack front for v/c. R = 0.2 
(delay time cir/h — 0.3177) 






5.9 Conclusions 


157 





K/faJxa) 


158 


Dynamic Analysis: Running Crack 



Figure 5.8: DSIF variation with time for different positions along the crack front for v/cr = 0.6 
(delay time c\rjh = 0.03177) 
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Figure 5.9: DSIF variation with time for different positions along the crack front for v/c R = 0.8 
(delay time cir/h = 0.03177) 
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Figure 5.10: DSIF variation along the crack front for different instants of time for v/c.j { = 0.4 
(delay time c\rjh — 0.032) 
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Figure 5.11: DSIF variation along the crack front for different instants of time for v/cr - 0.8 
(delay time c\r/h — 0.032) 
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Figure 5.12: DSIF variation along the crack front for different instants of time for v/cn - 0.8, 
using only DSIF values at corner nodes of the elements on crack front 
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Figure 5.14. DSIF variation with time at crack vertex for different speeds of crack propagation 
(delay time c\rjh = 0.3177) 
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Figure 5.15: DSIF variation with time at mid-surface for different speeds of crack propagation 
(delay time c\rjh = 0.3177) 




K r /(o 0 ^na) 


166 


Dynamic Analysis: Runni ng C rack 



Figure 5.16: DSIF variation with time at crack vertex for different speeds of crack propagation 
(delay time c lT /h = 0.03177) 
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Figure 5.17: DS1F variation with time at mid-surface for different speeds of crack propagation 
(delay time c\r/h = 0.03177) 
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Figure 5.21: Plots of Kj\t) with time for different positions along crack front for intersection 
angle of 101° (delay time c\r/h = 0.3177) 




Chapter 6 

Conclusions and Future Directions 


This chapter summarises this thesis, and will set forth goals for future work. 


6.1 Conclusions 

In this work, the effect of free surfaces on fracture characteristics for elastostatic and 
elastodynamic problems has been investigated using three-dimensional boundary element 
method. Free surface effect has been studied for a through center cracked finite body 
for three situations: quasi-static conditions, stationary crack under impact loading and 
running crack under impact loading. 

For finite geometries, finite element analysis is not able to obtain an accurate stress 
state adjacent to the crack tip in the region close to the free surface. The problems of 
brittle fracture involving intersection of the crack fronts with the free surfaces still requires 
extremely careful examination. The BEM has special suitability for problems of fracture 
mechanics. Its capability to accurately model steep stress gradients is the main motivation 
for its use in this study to analyse free surface effect for through cracks. 

Stress intensity factors from a boundary element analysis can be directly calculated 
from values of tractions at nodes on the crack front. Values of J" integral calculated along 
the crack front correspond well with the Kj values calculated from nodal traction values 
in quasi-static analysis. Thus, it is possible to relate the evaluated state of stress at the 
crack vertex, where crack front meets the free surface, with a value of SIF. This SIF value 
reflects the influence of the exponent of singularity existing at the crack vertex. A lower 
value of SIF is computed for a weaker singularity at the vertex and a higher value of SIF is 
computed for a stronger singularity at the vertex. For the same level of mesh refinement, 
value of SIF is seen to be sensitive to change in singularity of stresses at the vertex brought 
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about by changes in either Poisson’s ratio of the material or intersection angle of the 
crack front with the free surface. It becomes possible to estimate the intersection angle 
at which the SIF at the crack vertex is equal to or higher than mid-surface SIP value. So 
the intersection angle at which free surface vertex point becomes critical along with other 
parts of the crack front gets estimated. 

Values of DSIF calculated from nodal tractions at crack front, show sensitivity to loss of 
square-root dominance in transient fields. Lower values of DSIF an' found at. corner point 
for Mode-I loading for low speeds of crack propagation, confirming the existence of lower 
singularity in stresses at corner point. But at crack speeds beyond v — (). lc;,>, higher values 
of DSIF are computed compared to the rest of the crack front. The study correctly predicts 
that for high speeds of crack propagation, crack front maintains its straight profile. 

It should become possible to capture free surface ret ardation of crack growth accurately 
in different crack propagation problems with the methodology used in this paper for both 
static and dynamic conditions. 

6.2 Future Directions 

The presented methodology may be further tested on actual crack propagat ion problems. In 
this study, only self-similar crack propagation was considered. The presented methodology 
may be used to resolve the question of uniqueness of relationship between Kf and v in high 
speed crack propagation experiments and simulations. For dynamic problems, dynamic 
integral may be computed and its correspondence with DSIF values computed from nodal 
tractions on crack front verified. The work can be easily extended to study free surface 
effect for interfacial cracks. 
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